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Abstract 



An algorithm for detecting unstable periodic orbits in chaotic systems [Phys. Rev. 
E, 60 (1999), pp. 6172-6175] which combines the set of stabihsing transformations 
proposed by Schmelcher and Diakonos [Phys. Rev. Lett., 78 (1997), pp. 4733-4736] 
with a modified semi-implicit Euler iterative scheme and seeding with periodic orbits 
of neighbouring periods, has been shown to be highly efficient when applied to low- 
dimensional system. The difficulty in applying the algorithm to higher dimensional 
systems is mainly due to the fact that the number of stabilising transformations 
grows extremely fast with increasing system dimension. In this thesis, we construct 
stabilising transformations based on the knowledge of the stability matrices of al- 
ready detected periodic orbits (used as seeds). The advantage of our approach is 
in a substantial reduction of the number of transformations, which increases the 
efficiency of the detection algorithm, especially in the case of high-dimensional sys- 
tems. The dependence of the number of transformations on the dimensionality of 
the unstable manifold rather than on system size enables us to apply, for the first 
time, the method of stabilising transformations to high-dimensional systems. An- 
other important aspect of our treatment of high-dimensional flows is that we do not 
restrict to a Poincare surface of section. This is a particularly nice feature, since 
the correct placement of such a section in a high-dimensional phase space is a chal- 
lenging problem in itself. The performance of the new approach is illustrated by 
its application to the four-dimensional kicked double rotor map, a six- dimensional 
system of three coupled Henon maps and to the Kuramoto-Sivashinsky system in 
the weakly turbulent regime. 
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Chapter I 



Introduction 



The successes of the differential equation paradigm were impressive and 
extensive. Many problems, including basic and important ones, led to 
equations that could be solved. A process of self-selection set in, whereby 
equations that could not be solved were automatically of less interest 
than those that could. 
/. Stewart 

In this chapter we start in §1.1 by giving a brief primer into the theory of dynam- 
ical systems. Here our intention is not to give an exhaustive review (concise reviews 
on the subject are given in [33, 42, 68]). Rather, it is to illustrate the role played by 
periodic orbits in the development of the theory. In §1.2 a brief introduction to the 
periodic orbit theory is provided, followed by a discussion concerning the efficient 
detection of unstable periodic orbits (UPOs). Section 1.3 looks at the application 
to high-dimensional systems, in particular, large nonequilibrium systems that are 
extensively chaotic. It is well known that numerical methods can both introduce 
spurious chaos, as well as suppress it [10, 98]. Thus in §1.4 we discuss some of the 
numerical issues which can arise when detecting UPOs for a chaotic dynamical sys- 
tem. We give an overview of the objectives of this thesis in §1.5. The final section, 
1.6, details the contribution to the literature of this thesis. 

1.1 History, theory and applications 

Although the subject of modern dynamical systems has seen an explosion of interest 
in the past thirty years - mainly due to the advent of the digital computer - its roots 
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firmly belong at the foot of the twentieth century. Partly motivated by his work on 
the famous three body problem, the French mathematician and philosopher Henri 
Jules Poincare was to revolutionise the study of nonlinear differential equations. 

Since the birth of the calculus, differential equations have been studied both 
in their own right and for modeling phenomena in the natural sciences. Indeed, 
Newton considered them in his work on differential calculus [67]^ as early as 1671. 
One of the earliest examples of a first order equation considered by Newton was 

^ = 1 — 3x + y + x'^ + xy. (1.1) 
ax 

A solution of this equation for the initial condition y{0) = can be obtained as 
follows: start with 

y = + --- 

and insert this into Eq. (1.1); integrating yields 

y = X H , 

repeating the process with the new value of y gives 

y = X — x"^ + ■ ■ ■ . 

One can imagine continuing this process ad infinitum, leading to the following solu- 
tion of Eq. (1.1) 

2 "^4 "^"^ "^fi 

y = X — X H — X X H X x + ■ ■ ■ 

^ 3 6 30 45 

(for further details see [36]). 

The preceding example demonstrates one of the main differences between the 
classical study of differential equations and the current mindset. The classical study 
of nonlinear equations was local, in the sense that individual solutions where sought 
after. Most attempts in essence, involved either an approximate series solution or 
determining a transformation under which the equation was reduced either to a 
known function or to quadrature. 

In his work on celestial mechanics [74], Poincare developed many of the ideas 



^Newton's De Methodis Serierum et Fluxionum was written in 1671 but Newton failed to get 
it published and it did not appear in print until John Colson produced an English translation in 
1736. 
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underpinning modern dynamical systems. By working with sets of initial conditions 
rather than individual solutions, he was able to prove that complicated orbits existed 
for which no closed solution was available; Poincare had caught a glimpse of what 
is popularly coined "chaos" nowadays. 

Although there was continued interest from the mathematical community; most 
notably Birkhoff in the 1920s and the Soviet mathematicians in the 1940s - Kol- 
mogorov and students thereof - it was not until the 1960s that interest from the gen- 
eral scientific community was rekindled. In 1963 the meteorologist Edward N. Lorenz 
published his now famous paper "deterministic nonperiodic flow" [61] where a simple 
system describing cellular convection was shown to exhibit extremely complicated 
dynamics. Motion was bounded, displayed sensitivity to initial conditions and was 
aperiodic; Lorenz had witnessed the first example of a chaotic attractor. 

Around the same time, the mathematician Steve Smale was using methods from 
differential topology in order to prove the existence of a large class of dynamical 
systems (the so called axiom-A systems), which were both chaotic and structurally 
stable at the same time [91]. Along with examples such as the Lorenz model above, 
scientists where lead to look beyond equilibrium points and limit cycles in the study 
of dynamical processes. It became clear that far from being a mathematical oddity, 
the chaotic evolution displayed by many dynamical systems was of great practical 
importance. 

Today the study of chaotic evolution is widespread throughout the sciences where 
the tools of nonlinear analysis are used extensively. There remain many open ques- 
tions and the theory of dynamical systems has a bright and challenging future. The 
prediction and control [70, 86] of deterministic chaotic systems is an important area 
which has received a lot of attention over the past decade, whilst the extension of 
the theory to partial differential equations [79, 93] promises to give fresh insight into 
the modeling of fully developed turbulence. However, perhaps the most promising 
area of future research lies in the less mathematically minded disciplines such as 
biology, economics and the social sciences, to name a few. 

1.2 Periodic orbits 

Periodic orbits play an important role in the analysis of various types of dynam- 
ical systems. In systems with chaotic behaviour, unstable periodic orbits form a 
"skeleton" for chaotic trajectories [16]. A well regarded definition of chaos [23] re- 
quires the existence of an infinite number of UPOs that are dense in the chaotic 
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set. Different geometric and dynamical properties of chaotic sets, such as natural 
measure, Lyapunov exponents, fractal dimensions, entropies [69], can be determined 
from the location and stability properties of the embedded UPOs. Periodic orbits 
are central to the understanding of quantum-mechanical properties of nonseparable 
systems: the energy level density of such systems can be expressed in a semiclassical 
approximation as a sum over the UPOs of the corresponding classical system [35]. 
Topological description of a chaotic attractor also benefits from the knowledge of 
periodic orbits. For example, a large set of periodic orbits is highly constraining 
to the symbolic dynamics and can be used to extract the location of a generating 
partition [20, 73]. The significance of periodic orbits for the experimental study of 
dynamical systems has been demonstrated in a wide variety of systems [58], espe- 
cially for the purpose of controlling chaotic djTiamics [70] with possible application 
in communication [5]. 

1.2.1 Periodic orbit theory 

Briefly put, the periodic orbit theory provides a machinery which enables us to 
use the knowledge provided by the properties of individual solutions, such as their 
periods, location and stabilities, to make predictions about statistics, e.g. Lyapunov 
exponents, entropies, and so on. The dynamical systems to be discussed in this 
section are smooth n-dimensional maps of the form Xj+i = /(xj), where Xi is an 
n-dimensional vector in the ?7,-dimensional phase space of the system. 

Now, in order for the results to be quoted to hold, we assume that the attractor 
of / is both hyperbolic and mixing. A hyperbolic attractor is one for which the 
following two conditions hold: (i) there exist stable and unstable manifolds at each 
point of the attractor whose dimensions, and Till •) El'I'S the same for each point on 
the attractor, with + riu = n, and (ii) there exists a constant K > 1 such that 
for all points, x, on the attractor, if a vector u is chosen tangent to the unstable 
manifold, then 

\\Df{x)u\\>K\\u\\, (1.2) 
and if u is chosen tangent to the stable manifold 

\\Df{x)u\\<\\u\\/K. (1.3) 

Here Df{x) denotes the Jacobian matrix of the map / evaluated at the point x. By 



1.2 Periodic orbits 



5 



mixing we mean that for any two subsets Ai, A2 in the phase space, we have 

hm n f (A2)] = fi{A,)fi{A2), (1.4) 

j— >oo 

where fi is the natural measure of the attractor. In other words, the system will 
evolve over time so that any given open set in phase space will eventually overlap 
any other given region. 

Let us denote the magnitudes of the eigenvalues of the Jacobian matrix for the 
p times iterated map evaluated at the jth fixed point by Ai-,-, . . . , A„j. Suppose 
that the number of unstable eigenvalues, i.e. Xij > 1, is given by and further, 
that we order them as follows 

Alj > ■ ■ ■ > Xriuj > 1 > X{nu+l)j > ■ • • > Xnj- (1-5) 

Let Lj denote the product of unstable eigenvalues at the jth fixed point of 

Lj = \lj\2j ■ ■ ■ Xnuj- (1-6) 

Then the principal result of the periodic orbit theory is the following: given a subset 
A of phase space, one may define its natural measure to be 

/i(A) = lim fip{A), (1.7) 

where 

MA) = ^LT\ (1.8) 

j 

Here the sum is over all fixed points of in A; a derivation of Eq. (1.7) may be 
found in [34]. 

This result leads to several important consequences, for example, it can be shown 
that the Lyapunov numbers of / are given by 

log Xp = lim i V LJ^ log Xpj, (1.9) 
j 

whilst an analogous result exists for the topological entropy 

hr = lim -InA^p, (1.10) 

p^oo p 

where Np denotes the number of fixed points of the map These and similar 
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results obtained within the periodic orbit theory show that knowledge of the UPOs 
can yield a great deal of information concerning the properties of a chaotic dynamical 
system. Thus making their efficient detection highly desirable. For further details, 
a thorough review of the periodic orbit theory is given in the book by Cvitanovic et 
al [17]. 

At this stage, it is important to point out that most systems of interest turn 
out not to be hyperbolic, in particular, the dynamical systems studied in this thesis 
are non-hyperbolic. Hyperbolic systems, however, remain important due to the fact 
that they are more tractable from a mathematical perspective. Indeed, most rigorous 
results in dynamical systems are for the case of hyperbolic systems, and although 
much of the theory is believed to transfer over to the non-hyperbolic case there are 
very few rigorous results. 

1.2.2 Efficient detection of UPOs 

We have seen that the role of UPOs in chaotic systems is of fundamental theoretical 
and practical importance. It is thus not surprising that much effort has been put 
into the development of methods for locating periodic solutions in different types 
of dynamical systems. In a limited number of cases, this can be achieved due to 
the special structure of the systems. Examples include the Biham-Wenzel method 
applicable to Henon-like maps [3] , or systems with known and well ordered symbolic 
dynamics [41]. For generic systems, however, most methods described in the liter- 
ature use some type of an iterative scheme that, given an initial condition (seed), 
converges to a periodic orbit of the chaotic system. In order to locate all UPOs with 
a given period p , the convergence basin of each orbit for the chosen iterative scheme 
must contain at least one seed. The seeds are often chosen either at random from 
within the region of interest, from a regular grid, or from a chaotic trajectory with 
or without close recurrences. Typically, the iterative scheme is chosen from one of 
the "globally" convergent methods of quasi-Newton or secant type. However, expe- 
rience suggests that even the most sophisticated methods of this type suffer from 
a common problem: with increasing period, the basin size of the UPOs becomes 
so small that placing a seed within the basin with one of the above listed seeding 
schemes is practically impossible [64]. 

A different approach, which appears to effectively deal with the problem of re- 
duced basin sizes has been proposed by Schmelcher and Diakonos (SD) [83, 84] . 
The basic idea is to transform the dynamical system in such a way that the UPOs 
of the original system become stable and can be located by simply following the 
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evolution of the transformed dynamical system. That is, to locate period-p orbits 
of a discrete dynamical system 

U: x,+i = /(x,), /:R"^M\ (1.11) 

one considers an associated flow 

E: j-^=Cg{x), (1.12) 

where g{x) = f^{x) — x and C is an x n constant orthogonal matrix. It is easy 
to see that the map f^{x) and flow E have identical sets of fixed points for any C, 
while C can be chosen such that unstable period-p orbits of U become stable fixed 
points of E. 

Since it is not generally possible to choose a single matrix C that would stabilise 
all UPOs of U, the idea is to find the smallest possible set of matrices C = {Cfc}^^, 
such that, for each UPO of U, there is at least one matrix C G C that transforms 
the unstable orbit of U into a stable fixed point of E. To this end, Schmelcher and 
Diakonos have put forward the following conjecture [83] 

Conjecture 1.1. Let Csd be the set of all n x n orthogonal matrices with only ±1 
non-zero entries. Then, for any n x n non-singular real matrix G, there exists a 
matrix C G Csd such that all eigenvalues of the product CG have negative real parts. 

Observation 1.1. The set Csd forms a group isomorphic to the Weyl group Bn [48], 
i.e. the symmetry group of an n-dimensional hypercube. The number of matrices in 
Csd is K = 2''n\. 

The above conjecture has been verified for n < 2 [85], and appears to be true for 
n > 2, but, thus far, no proof has been presented. According to this conjecture, any 
periodic orbit, whose stability matrix does not have eigenvalues equal to one, can 
be transformed into a stable fixed point of E with G G Csd- In practice, to locate 
periodic orbits of the map U, we try to integrate the flow E from a given initial 
condition (seed) using different matrices from the set Csd- Some of the resulting 
trajectories will converge to fixed points, while others will fail to do so, either leaving 
the region of interest or failing to converge within a specified number of steps. 

The main advantage of the SD approach is that the convergence basins of the 
stabilised UPOs appear to be much larger than the basins produced by other iter- 
ative schemes [21, 55, 84], making it much easier to select a useful seed. Moreover, 
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depending on the choice of the stabihsing transformation, the SD method may con- 
verge to several different UPOs from the same seed. 

The flow S can be integrated by any off-the-shelf numerical integrator. Schmelcher 
and Diakonos have enjoyed considerable success using a simple Euler method. How- 
ever, the choice of integrator for this problem is governed by considerations very 
different from those typically used to construct an ODE solver. Indeed, to locate 
a fixed point of the flow, it may not be very efficient to follow the flow with some 
prescribed accuracy. Therefore, local error considerations, for example, are not as 
important. Instead, the goal is to have a solver that can reach the fixed point in as 
few integration steps as possible. In fact, as shown by Davidchack and Lai [19], the 
efficiency of the method can be improved dramatically when the solver is constructed 
specifically with the above goal in mind. In particular, recognizing the typical stiff- 
ness of the flow S, Davidchack and Lai have proposed a modified semi-implicit Euler 
method 

= x, + [l3s,C'^ - G.]-^g{x,) , (1.13) 

where /? > is a scalar parameter, Si = \ \g{xi)\\ is an L2 norm, Gi = Dg{xi) is the 
Jacobian matrix, and "T" denotes transpose. Note that, away from the root of the 
above iterative scheme is a semi-implicit Euler method with step size h = 
and, therefore, can follow the flow S with a much larger step size than an explicit 
integrator (e.g. Euler or Runge-Kutta). Close to the root, the proposed scheme can 
be shown to converge quadratically [55] , analogous to the Newton- Raphson method. 

Another important ingredient of the algorithm presented in [19] is the seed- 
ing with already detected periodic orbits of neighbouring periods. This seeding 
scheme appears to be superior to the typically employed schemes and enables fast 
detection of plausibly all^ periodic orbits of increasingly larger periods in generic 
low-dimensional chaotic systems. For example, for the Ikeda map at traditional 
parameter values, the algorithm presented in [19] was able to locate plausibly all 
periodic orbits up to period 22 for a total of over 10^ orbit points. Obtaining a com- 
parable result with generally employed techniques requires an estimated 10^ larger 
computational effort. 

While the stabilisation approach is straightforward for relatively low- dimensional 
systems, direct application to higher-dimensional systems is much less efficient due 
to the rapid growth of the number of matrices in Csd- Even though it appears that, 
in practice, far fewer transformations are required to stabilise all periodic orbits of 
a given chaotic system [72], the sufficient subset of transformations is not known a 

2See §1.4 
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priori. It is also clear that the route of constructing a universal set of transformations 
is unlikely to yield substantial reduction in the number of such transformations. 
Therefore, a more promising way of using stabilising transformations for locating 
periodic orbits in high-dimensional systems is to design such transformations based 
on the information about the properties of the system under investigation. 

1.3 Extended systems 

The periodic orbit theory is well developed for low-dimensional chaotic dynamics - at 
least for axiom-A systems [17]. The question naturally arises as to whether or not the 
theory has anything to say for extended systems. At first glance the transition from 
low-dimensional chaotic dynamics to fully developed spatiotemporal chaos may seem 
rather optimistic. However, recent results have shown that certain classes of PDEs 
turn out to be less complicated than they initially appear, when approached from a 
dynamical systems perspective. Indeed, under certain conditions their asymptotic 
evolution can be shown to lie on a finite dimensional global attractor [78, 79, 93]. 
Further, by restricting to equations of the form 



where A is a linear differential operator, an even stronger result may be obtained. 
Such equations are termed evolution equations and their asymptotic dynamics can 
be shown to lie on a smooth, finite dimensional manifold, known as the inertial 
manifold [78]. In contrast to the aforementioned global attractor which may have 
fractal like properties this leads to a complete description of the dynamics by a finite 
number of modes; higher modes being contained in the geometrical constraints which 
define the manifold. 

A variety of methods for determining all UPOs up to a given length exist for low- 
dimensional dynamical systems (see Chapter 2). For more complex dynamics, such 
as models of turbulence in fluids, chemical reactions, or morphogenesis in biology 
with high - possibly infinite - dimensional phase spaces, such methods quickly run 
into difficulties. The most computationally demanding calculation to date, has been 
performed by Kawahara and Kida [52]. They have reported the detection of two 
three-dimensional periodic solutions of turbulent plane Couette flow using a 15, 422- 
dimensional discretisation, whilst more recently Viswanath [95] has been able to 
detect both periodic and relative periodic motions in the same system. It is hoped 




(1.14) 
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that such solutions may act as a basis to infer the manner in which transitions to 
turbulence can occur. 

Our goal is somewhat more modest. We will apply our method to the model 
example of an extended system which exhibits spatiotemporal chaos; the Kuramoto- 
Sivashinsky equation 



It was first studied in the context of react ion- diffusion equations by Kuramoto and 
Tsuzuki [56], whilst Sivashinsky derived it independently as a model for thermal 
instabilities in laminar flame fronts [90]. It is one of the simplest PDEs to exhibit 
chaos and has played a leading role in studies on the connection between ODEs and 
PDEs [27, 49, 54]. 

It is the archetypal equation for testing a numerical method for computing 
periodic solutions in extended systems, and has been considered in this context 
in [9, 57, 100], where many UPOs have been detected and several dynamical aver- 
ages computed using the periodic orbit theory. Note that the attractor of the system 
studied in [9] is low dimensional, whilst those studied in [57, 100] have higher in- 
trinsic dimension. Recently the closely related complex Ginzburg-Landau equation 



has been studied within a similar framework [60], where Eq. (1.16) is transformed 
into a set of algebraic equations which are then solved using the Levenberg-Marquardt 
algorithm. 

1.4 A note on numerics 

Often the physical models put forward by the applied scientists are extremely com- 
plex and, thus, not open to attack via analytical methods. This necessitates the use 
of numerical simulations in order to analyse and understand the models - particu- 
larly in the case where chaotic behaviour is allowed. However, in those cases one 
can always wonder what one is really computing, given the limitations of floating 
point systems. This leads to the important question of whether or not the computed 
solution is "close" to a true solution of the system of interest. In the case of locating 
a periodic orbit on a computer, we would like to know whether the detected orbit 
actually exists in the real system. It is well known that in any numerical calcula- 




(1.15) 




(1.16) 
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tion accuracy is limited by errors due to roundoff, discretisation and uncertainty of 
input data [92]. Tlie difficulty here, lies in the fact that the solutions of a chaotic 
dynamical system display extreme sensitivity upon initial conditions, thus, any tiny 
error will result in the exponential divergence of the computed solution from the 
true one. 

For discrete hyperbolic systems, an answer to the question of validity is provided 
by the following shadowing lemma due to Anosov and Bowen [1, 7] 

Lemma 1.1. Given a discrete hyperbolic system, 

Xi+i = /(xi), (1.17) 

then Ve > 0, 3(5 > such that every 6 -pseudo-orbit for f is e-shadowed by a unique 
real orbit. 

By 5-pseudo-orbit, we mean a computed sequence (pi)igz such that 

\pi+,-f{pi)\<6, (1.18) 

that is to say, roundoff error at each step of the numerical orbit is bounded above 
by 6. Such an orbit is said to be e-shadowed if there exists a true orbit (xj)jgz such 
that 

\xi-pi\<e VieZ. (1.19) 

Unfortunately the class of hyperbolic systems is highly restrictive since such 
systems are rarely encountered in real problems. For non-hyperbolic systems - such 
as those studied in this work - shadowing results are limited to low dimensional 
maps [11, 39]. Even then, shadowing can only be guaranteed for a finite number of 
steps N, which is likely to be a function of the system parameters. Further, it can 
be shown that trajectories of non-hyperbolic chaotic systems fail to have long time 
shadowing trajectories at all, when unstable dimension variability persists [22, 24, 
81, 82]. Although the idea of shadowing goes some way towards making sense of 
numerical simulations of chaotic systems, it does not answer the question of whether 
the numerical orbit corresponds to any real one. Therefore, we need tools which can 
rigorously verify the existence of the corresponding real periodic orbits. 

Such methods may also be used to determine the completeness of sets of periodic 
orbits, however, this requires that the entire search is conducted using rigorous 
numerics and this approach is inefficient for generic dynamical systems. In general, 
it is not possible to prove, within our approach, the completeness of the detected 
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sets of UPOs. Rather, a stopping criteria must be deduced after which we can 
say, with some certainty, that all UPOs of period-p have been found; our assertion 
of completeness will be based upon the plausibility argument. The following three 
criteria are used for the validation of the argument 

(i) Methods based on rigorous numerics (e.g. in [28]) have located the same UPOs 
in cases where such comparison is possible (usually for low periods, since these 
methods are less efficient). 

(ii) Our search strategy scales with the period p (see §3.4 and [21]). If we can tune 
it to locate all UPOs for low periods (where we can verify the completeness 
using (i)), it is hkely (but not provably) capable of locating all UPOs of higher 
periods as well. 

(iii) For maps with symmetries, we test the completeness by verifying that all the 
symmetric partners for all detected UPOs have been found (see §3.4.1 and 
§3.4.2). 

Of course, the preceding discussion can only be applied to discrete systems. 
1.4.1 Interval arithmetic 

There are several approaches towards a rigorous computer assisted study of the 
existence of periodic orbits. Most make use of the Brouwer fixed point theorem [8], 
which states that if a convex, compact set X C M" is mapped by a continuous 
function / into itself then / has a fixed point in X. Such rigorous methods tend to 
fall into one of two classes: (i) topological methods based upon the index properties 
of a periodic orbit or (ii) interval methods. In our discussion we restrict attention 
to interval methods since they are the most common in practice. Techniques based 
on the index properties are in general less efficient; although recent work has seen 
the ideas extended to include infinite dimensional dynamical systems [30, 99], in 
particular, in [99] several steady states for the Kuramoto-Sivashinsky equation have 
been verified rigorously. 

Interval methods are based on so called interval arithmetic - an arithmetic de- 
fined on sets of intervals [65]. Any computation carried out in interval arithmetic 
returns an interval which is guaranteed to contain both the true solution and the 
numerical one. Thus, by using properly rounded interval arithmetic, it is possible to 
obtain rigorous bounds on any numerical calculations. In what follows an interval 
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is defined to be a compact set x C M, i.e. 

X = [a, 6] = {x : a < X < 6}, 

where we use boldface letters to denote interval quantities and lowercase maths 
italic to denote real quantities. By an n-dimensional interval vector, we refer to the 
ordered n-tuple of intervals v = {xi, . . . ,x„}. Note that this leads readily to the 
definition of higher dimensional objects. 

Arithmetic on the set of intervals is naturally defined in the following way: let 
us denote by o one of the standard arithmetic operations +, — , ■ and /, then the 
extension to arbitrary intervals Xi and X2 must satisfy the condition 

Xl O X2 = {X = Xl O X2 : Xl e Xl, X2 G X2}, 

where, in the case of division, the interval X2 must not contain the number zero. 
Importantly, the resulting interval is always computable in terms of the endpoints, 
for example, let xi = [a, b] and X2 = [c, d] then the four basic arithmetic operations 
are given by 



Xl + X2 = 


[a + c,b + d], 


(1.20) 


Xl - X2 = 


[a — d,b — c], 


(1.21) 


Xl * X2 = 


[minjac, ad, be, bd}, max{ac, ad, be, bd}]. 


(1.22) 


1/xi = 


[1/6, 1/a], 


(1.23) 


X1/X2 = 


Xl * I/X2. 


(1.24) 



This allows one to obtain bounds on the ranges of real valued functions by writing 
them as the composition of elementary operations. For example, if 

/(x) = x(x - 1), 

then 

/([0,1]) = [0,1]([0,1]-1) = [-1,0], 

note the exact range [—1/4,0] C [—1,0] as expected. 

Combined with the Brouwer fixed point theorem, interval arithmetic enables 
us to prove the existence of solutions to nonlinear systems of equations. In §1.2.2 
we saw that the periodic orbit condition is equivalent to the following system of 
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nonlinear equations 



g{x) = 



where g{x) = f^{x) — x. In order to investigate the zeros of the function g one may 
apply the Newton operator to the n-dimensional interval vector x 



where g'{'x) is the interval matrix containing all Jacobian matrices of the form g'{x) 
for x G X, and Xq is an arbitrary point belonging to the interval x. Applying the 
Brouwer fixed point theorem in the context of the Newton interval operator leads 
to the following Theorem. 

Theorem 1.1. // N{x) C int(a;) then g{x) = has a unique solution in x. If 
N{x) n cc = then there are no zeros of g in x. 

For a proof, see for example, [28]. 

In practice, the following algorithm may be applied to verify the existence of a 
numerical orbit: (i) start by surrounding the orbit by an n-dimensional interval of 
width e, where e is an integer multiple of the precision, e, with which the orbit is 
known, (ii) then apply the Newton operator to the interval, if A^(x) C x there is 
exactly one orbit in x, else if iV(x) fix = no orbit of g lies in x, (iii) if neither of the 
above hold then either the orbit is not a true one, or else, e needs to be increased. 

In [28, 29] interval arithmetic has been applied to various two-dimensional maps, 
note however, that in applications the Newton operator is replaced by the following 
method due to Krawczyk 



here A is a preconditioning matrix. The Krawczyk operator of Eq. (1-26) has the 
advantage that it does not need to compute the inverse of g', thus it can be used for 
a wider class of systems than the Newton operator. 

1.5 Overview 

In this thesis, we present an extension of the method of stabilising transformations 
to high- dimensional systems. Using periodic orbits as seeds, we construct stabilising 
transformations based upon our knowledge of the respective stability matrices. The 
major advantage of this approach as compared with the method of Schmelcher and 



A^(x) = xo - {g'{^)) ^9{xo), 



(1.25) 



K{x.) =xo- Ag{xo) - {Ag'ip^) - J)(x - xq). 



(1.26) 
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Diakonos is in a substantial reduction of the number of transformations. Since in 
practice, high-dimensional systems studied in dynamical systems typically consist of 
low-dimensional chaotic dynamics embedded within a high- dimensional phase space, 
we are able to greatly increase the efficiency of the algorithm by restricting the 
construction of transformations to the low- dimensional dynamics. An important 
aspect of our treatment of high-dimensional flows is that we do not restrict to a 
Poincare surface of section (PSS). This is a particularly nice feature, since the phase 
space topology for a high-dimensional flow is extremely complex, and the correct 
placement of such a surface is a nontrivial task. 

In Chapter 2 we review common techniques for detecting UPOs, keeping with 
the theme of the present work our arrangement is biased towards those methods 
which are readily applicable in higher dimensions. We begin Chapter 3 by introduc- 
ing the method of stabilisation transformations (ST) in its original form. In §3.2 
we study the properties of the STs for n = 2. We extend our analysis to higher 
dimensional systems in §3.3, and show how to construct STs using the knowledge 
of the stabihty matrices of already detected periodic orbit points. In particular, we 
argue that the stabilising transformations depend essentially on the signs of unsta- 
ble eigenvalues and the directions of the corresponding eigenvectors of the stability 
matrices. Section 3.4 illustrates the application of the new STs to the detection of 
periodic orbits in a four- dimensional kicked double rotor map and a six- dimensional 
system of three coupled Henon maps. In Chapter 4 we propose and implement an 
extension of the method of STs for detecting UPOs of flows as well as unstable 
spatiotemporal periodic solutions of extended systems. We will see that for high- 
dimensional flows - where the choice of PSS is nontrivial - it will pay to work in 
the full phase space. In §4.1 we adopt the approach often taken in subspace itera- 
tion methods [62], we construct a decomposition of the tangent space into unstable 
and stable orthogonal subspaces, and construct STs without the knowledge of the 
UPOs. This is particularly useful since in high dimensional systems it may prove 
difficult to detect even a single periodic orbit. In particular, we show that the use 
of singular value decomposition to approximate the appropriate subspaces is prefer- 
able to that of Schur decomposition, which is usually employed within the subspace 
iteration approach. The proposed method for detecting UPOs is tested on a large 
system of ODEs representing odd solutions of the Kuramoto-Sivashinsky equation 
in §4.2. Chapter 5 summarises this work and looks at further work that should be 
undertaken to apply the methods presented to a wider range of problems. 
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1.6 Thesis results 

The main results of this thesis are pubhshed in [13, 14, 15], the important points of 
which are detailed below. 

1) Efficient detection of periodic orbits in chaotic systems by stabilising 
transformations. 

• A proof of Conjecture 1.1 for the case = 2 is presented. In other words, 
we show that any two by two matrix may be stabilised by at least one matrix 
belonging to the set proposed by Schmelcher and Diakonos. 

• Analysis of the stability matrices for the two-dimensional case is provided. 

• The above analysis is used to construct a smaller set of stabilising transfor- 
mations. This enables us to efficiently apply the method to high-dimensional 
systems. 

• Experimental evidence is provided showing the successful application of the 
new set of transformations to high- dimensional {n > 4) discrete dynamical 
systems. For the first time, plausibly complete sets of periodic orbits are 
detected for high-dimensional systems. 

2) On the use of stabilising transformations for detecting unstable peri- 
odic orbits for the Kuramoto-Sivashinsky equation. 

• The extension of the method of stabilising transformations to large systems of 
ODEs is presented. 

• We construct stabilising transformations using the local stretching factors of an 
arbitrary - not periodic - point in phase space. This is particularly important, 
since for very high- dimensional systems, finding small sets of UPOs to initiate 
the search becomes increasingly difficult. 

• The number of such transformations is shown to be determined by the system's 
dynamics. This contrasts to the transformations introduced by Schmelcher and 
Diakonos which grow with system size. 

• In contrast to traditional methods we do not use a Poincare surface of section, 
rather, we supply an extra equation in order to determine the period. 
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• Experimental evidence for the applicability of the aforementioned scheme is 
provided. In particular, we are able to calculate many time-periodic solutions 
of the Kuramoto-Sivashinsky equation using a fraction of the computational 
effort of generally employed techniques. 



Chapter II 



Conventional techniques for 
detecting periodic orbits 

Science is built up of facts, as a house is with stones. But a collection of 
facts is no more a science than a heap of stones is a house. 
H. J. Poincare 

The importance of efficient numerical schemes to detect periodic orbits has been 
discussed in the Introduction, where we have seen that the periodic orbits play 
an important role in our ability to understand a given dynamical system. In the 
following chapter we give a brief review of the most common techniques currently 
in use. In developing numerical schemes to detect unstable periodic orbits (UPO) 
there is much freedom. Essentially, the idea is to transform the system of interest 
to a new dynamical system which possesses the sought after orbit as an attracting 
fixed point. Most methods in the literature are designed to detect UPOs of discrete 
systems, the application to the continuous setting is then made by the correct choice 
of Poincare surface of section (PSS). For that reason in this chapter, unless stated 
otherwise, the term dynamical system will refer to a discrete dynamical system. 

2.1 Special cases 

In a select number of cases, efficient methods may be designed based on the special 
structure inherent within a particular system. In this section we discuss such meth- 
ods, with particular interest in the method due to Biham and Wenzel [3] applicable 
to Henon-type maps. In Chapter 3 we apply our method to a system of coupled 
Henon maps and validate our results against a method which is an extension of the 
Biham-Wenzel method. 
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2.1.1 One-dimensional maps 

Perhaps one of the simplest methods to detect UPOs in one-dimensional maps is 
that of inverse iteration. By observing that the unstable orbits of a one-dimensional 
map are attracting orbits of the inverse map, one may simply iterate the inverse map 
forward in time in order to detect UPOs. Since the inverse map is not one-to-one, at 
each iteration we have a choice of branch to make. By choosing the branch according 
to the symbolic code of the orbit we wish to find, we automatically converge to the 
desired cycle. 

The method cannot be directly applied to higher-dimensional systems since they 
typically have both expanding and contracting directions. However, if in the con- 
tracting direction the chaotic attractor is thin enough so as to be treated approxi- 
mately as a zero-dimensional object, then it may be possible to build an expanding 
one-dimensional map by projecting the original map onto the unstable manifold and 
applying inverse iteration to the model system. Orbits determined in this way will 
typically be "close" to orbits of the full system, and may be used to initiate a search 
of the full system using more sophisticated routines. 

Methods can also be constructed due to the fact that for one-dimensional maps 
well ordered symbolic dynamics exists. For ease of exposition, we shall describe one 
such method in the case of a unimodal mapping, /, that is, a mapping of the unit 
interval such that /(O) = /(I) = 0, f'{c) = and /"(c) < 0, where c G (0, 1) is the 
unique turning point of /. 

The symbolic dynamics description for a point x e [0, 1] is given by {sk} where 

s,^l^ if f^'-'K^) > ^c, ^2.1) 
1 otherwise. 

Here Xc is the unique turning point of the map /. Note that the order along the 
X-axis of two points x and y can be determined from their respective itineraries {sk} 
and {s'f,}. To see this, let us define the well ordered symbolic future 7 of the point x 
to be 

00 

7(,S) = 0.TO••• = 5]^l2-^ (2.2) 

k=l 

where S denotes the symbolic code of the point x and 

k 

= (mod2). 



2.1 Special cases 



20 



Now suppose that Si = s[, . . . , Sk = s'j^ and s^+i = and s'^_,_]^ = 1. Then it can be 
shown that 

fc+i 

^ Si (mod2) = 0. (2.3) 

i=l 

For a proof see for example [18]. 

Thus in order to detect UPOs of the map / one begins by determining the 
symbolic value 7c for the orbit S = ■ ■ ■ S1S2 ■ ■ ■ SpSiS2 ■ ■ ■ Sp - ■ ■ . Choosing a starting 
point Xq with symbolic value 7 from the unit interval, one may update the starting 
point by comparing its symbolic value, 7, against 7c, the value for the cycle. Using 
a binary search, this procedure will quickly converge to the desired orbit. 

The method can also be extended to deal with certain two-dimensional systems, 
in that case one must also define the well ordered symbolic past in order to uniquely 
identify orbits of the system. This idea has been applied to a number of different 
models, such as the Henon map, different types of billiard systems and the diamag- 
netic Kepler problem, to name a few [41]. 

We conclude this section by mentioning that for one- dimensional maps it is 
always possible to determine UPOs as the roots to the nonlinear equation g = 
f^{x) — X. Since it is straightforward to bracket the roots of a nonlinear equation in 
one dimension and thus apply any of a number of solvers to detect UPOs. We shall 
discuss methods for solving nonlinear equations in some detail in §2.2. 

2.1.2 The Biham-Wenzel method 

The Biham-Wenzel (BW) method has been developed and successfully applied in the 
detection of UPOs for Henon- type systems [3, 4, 75]. It is based on the observation 
that for maps such as the Henon map there exists a one-to-one correspondence 
between orbits of the map and the extremum configurations of a local potential 
function. For ease of exposition, we describe the method in the case of the Henon 
map which has the following form 

Xj+i = 1 — ax'^ + (2.4) 

when expressed as a one-dimensional recurrence relation. 

In order to detect a closed orbit of length p for the Henon map, we introduce a 
p-dimensional vector field, v{x), which vanishes on the periodic orbit 



— = Vi{x) = Xj+i — 1 -|- ax'^ — bxi-i, i = 1, . . . ,p. 



(2.5) 
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Now for fixed Xj+i, the equation Vi{x) = has two solutions which may be 
viewed as representing extremal points of a local potential function 

= i^Viix), Vi{x) = Xi{xi+i - hxi-i - 1) + -x]. (2.6) 

OXi 6 

Assuming the two extremal points to be real, one is a local minimum of Vi{x) and 
the other is a local maximum. The idea of BW was to integrate the flow (2.5) with 
an essential modification of the signs of its components 

— = SiVi, 1 = 1, ...,p, (2.7) 

where Sj = ±1. Note that Eq. (2.7) is solved subject to the periodic boundary 
condition Xp+i = xq. 

Loosely speaking, the modified flow will be in the direction of the local maximum 
of Vi{x) if Si = +1, or in the direction of the local minimum if Sj = —1. The 
differential equations (2.7) then drive an approximate initial guess towards a steady 
state of (2.5). Since the potential defined in Eq. (2.6) is unbounded for large \xi\, 
the flow will diverge for initial guesses far from the true trajectory. However, the 
basins of attraction for the method are relatively large, and it can be shown that 
convergence is achieved for all initial conditions as long as I ^2 1, 1 — 1, . . . are 
small with respect to i/a. For the standard parameter values a = 1.4, b = 0.3, BW 
report the detection of all UPOs for p < 28. 

An additional feature of the BW method is that the different sequences, {si}, 
when read as a binary code, turn out to be related to the symbolic code of the UPOs. 
Consider a periodic configuration {x,}, and the corresponding sequence {sj}. Then 
if we define 

Si = {-iysi, i = l,...,p, (2.8) 

it can be seen that for most trajectories, the sequence Si coincides with the symbolic 
dynamics Si of the Henon map, which we define as 5^ = 1 if > and Si = — 1 
if Xi < 0. Further, for a particular sequence {sj} the corresponding UPO does not 
always exist. In this case, the solution of (2.7) diverges, thus in principle the BW 
method detects all UPOs that exist and indicates which ones do not. 

To conclude, the BW method is an efficient method for detecting UPOs of Henon- 
type maps. However, it has been shown that for certain parameter values the method 
may fail to converge in some cases [31, 40]. In [31] it was shown that failure occurs 
in one of two ways: either the solution of Eq. (2.7) converges to a limit cycle rather 
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than a steady state, or uniqueness fails, i.e. two sequences converge to the same 
UPO. Extensions of the method include detection of all 2^ UPOs for the complex 
Henon map [4], as well as detection of UPOs for systems of weakly coupled Henon 
maps [75]. 

2.2 UPOs: determining roots of nonlinear func- 
tions 

One of the most popular detection strategies is to recast the search for UPOs in 
phase space into the equivalent problem of detecting zeros of a highly nonlinear 
function 

g(x) = fP{x) - X, X e R'^, (2.9) 

where p G Z''" denotes the period. This enables us to use tools developed for root 
finding to aid in our search. These methods usually come in the form of an iterative 
scheme 

Xi+i = $(xi), ^ = 0,1,2,..., (2.10) 

which, under certain conditions and for a sufficiently good guess, xq, guarantee 
convergence [92] . Note that the concern of the current work is with high- dimensional 
systems, thus we will not discuss those methods which are not readily applicable in 
this case. In particular we do not deal with those approaches designed primarily for 
one-dimensional systems. Examples and further references concerning root finding 
for one-dimensional systems may be found, for example, in [77]. 

2.2.1 Bisection 

Perhaps the most basic technique for detecting zeros of a function of one variable 
is the bisection method. Let : M — M be continuous on the interval [a,b], then 
if g{a) < and g{b) > 0, it is an immediate consequence of the intermediate value 
theorem that at least one root of g must lie within the interval (a, b). Assuming the 
existence of a unique zero of g & (a, b), we may proceed as follows: set Xi = |(a + 6) 
and evaluate g{xi), if g{xi) > then the root lies in the interval (a, Xi), otherwise it 
lies in the interval {xi,b). The iteration of this process leads to a robust algorithm 
which converges linearly in the presence of an isolated root. Problems may however 
occur when - as in our case - the detection of several roots is necessary. 

An extension of the bisection method to higher dimensional problems known 
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Figure 2.1: (a) The polyhedron ABDC is non-characteristic while the polyhedron 
AEDC is characteristic, (b) Application of the CB method to AEDC leads to suc- 
cesive characteristic polyhedra GEDC and HEDC. 



as characteristic bisection (CB) has recently been put forward by Vrahatis [6, 96]. 
Before giving a description of the CB method, we need to define the concept of a 
characteristic polygon. Let us denote by S the set of all n-tuples consisting of ±l's. 
Clearly, the number of distinct elements of S is given by 2". Now we form the 2" x n 
matrix A„ so that its rows are the elements of 5* without repetition. Next, consider 
an oriented ra-polyhedron, n„, with vertices Xi,i = 1,2, ... ,2'"'. We may construct 
the associated matrix S{g; Hn), whose ith row is given by 

sgng{xi) = (sgngiixi), sgng2ixi), sgngnixi)f , (2.11) 

where sgn denotes the sign function and g is given by Eq. (2.9). The polyhedron 
Un is said to be a characteristic polyhedron if the matrix S{g; n„) is equal to A„ 
up to a permutation of rows. A constructive process for determining whether a 
polyhedron is characteristic, is to check that all combination of signs are present at 
its vertices. In Figure 2.1a, one can immediately conclude that the polygon ABDC 
is not characteristic, whilst the polygon AEDC is. 

In our discussion we assume that we know the whereabouts of an isolated root 
of Eq. (2.9); in practice rigorous techniques such as those discussed in §1.4 can be 
used in order to locate the so called inclusion regions to initiate the search. The 
idea of the CB method is to surround the region in phase space which contains the 
root by a succession of n-polyhedra, u!h\ At each stage the polyhedra are bisected 



2.2 UPOs: determining roots of nonlinear functions 



24 



in such a way that the new polyhedra maintain the quahty of being characteristic. 
This refinement process is continued until the required accuracy is achieved. The 
idea is illustrated for the two-dimensional case in Figure 2.1b, where three iterates 
of the process are shown: = AEDC, H^^^ = GEDC and H^^^ = HEDC. 

Similar to the one-dimensional bisection, CB gives an unbeatably robust algo- 
rithm once an inclusion region for the UPO has been determined. It is independent 
of the stability properties of the UPO and guarantees success as long as the initial 
polyhedron is characteristic. However, there lies the crux of the method. In general 
to initiate the search, one needs to embed 2" points into an n-dimensional phase 
space, such that all combination of ±l's are represented. This is clearly a nontrivial 
task for high-dimensional problems such as those studied in this work. 

2.2.2 Newton- type methods 

In this section we discuss a class of iterative schemes which are popular in practice 
- the Newton- Raphson (NR) method and variants thereof. Due to excellent conver- 
gence properties, namely that convergence is ensured for sufficiently good guesses 
and that in the linear neighbourhood quadratic convergence is guaranteed, NR is 
the method of choice for many applications. 

Newton-Raphson method 

Let X* be a solution of the equation g{x*) = P{x*) — x* = 0, and assume that we 
have an initial guess Xi sufficiently close to x*, i.e. x* — Xi = 6xi where 6xi is small. 
Replacing x* in the above equation gives 



which, after Taylor series expansion of / and some algebraic manipulation yields 



where Df^{xi) is the Jacobian of f^{x) evaluated at x = Xj. Restricting to the linear 
part of Eq. 2.13 and setting Xj+i = Xi + 6xi, leads to the familiar Newton-Raphson 
algorithm 




(2.12) 



5X, + (Drix,) - ln)-\nx^) - + 0(5X^) = 0, 



(2.13) 



Xj+i = X,- {DfP{Xi) - In) \F{Xi) - Xi), 




(2.14) 
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In the case n = 1, Newton's method has the following geometrical interpretation: 
given an initial point Xj, we approximate g by the linear function 

L{x) = g{x,;) + Dg{xi){x - Xi), (2.15) 

which is tangent to g{x) at the point xi. We then obtain the updated guess Xj+i by 
solving the linear Eq. (2.15); see Figure 2.2. 

In order to speed up convergence, Householder proposed the following iterative 
scheme [47] 

Xi+, =x, + {k + l) (^Y/^y^j , A; = 1, 2, 3, ... . (2.16) 

Note that this is a generalization of the Newton algorithm applicable in the one- 
dimensional setting, where one keeps k + 1 terms in the Taylor series expansion of 
/. For k = 0, Newton's method is restored, whilst k = 1, gives the following third 
order method first discovered by the astronomer E. Halley in 1694 

Xi+i = Xi — f'^ -. (2-17) 
2(^0 - 99" 

The use of NR to detect periodic orbits is limited mainly due to the unpredictable 
nature of its global convergence properties, in particular for high dimensions. Typ- 
ically the basins of attraction will not be simply connected regions, rather, their 
geometries will be given by complicated, fractal sets, rendering any systematic de- 
tection strategy useless. Nevertheless, the method is unbeatably efficient when a 
good initial guess is known, or for speeding convergence in inaccurately determined 
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zeros. 



Multiple shooting algorithms 



Multiple shooting is a variant of the NR method, and has been developed for de- 
tecting period-p orbits for maps [17], particularly for increasing period, where the 
nonlinear function g becomes increasingly complex and can fluctuate excessively. In 
those cases the periodic orbit is more easily found as a zero of the following n x p 
dimensional vector function 



F(x) 



/(Xi) - X2 
fix2) - X3 



(2.18) 



f{X'p) X\ 

We now write Newton's method in the following more convenient form 

DF{x.)5x = -F(x), 



(2.19) 



where we define the Newton step 5x = Xj+i — Xj, and DF{x.) is the np x np matrix 



fix,) 

f{x2) -Ir 



f'{Xp_i) -1^ 



fix,) J 



(2.20) 



where I„ is the n x n identity matrix. Note that, due to the sparse nature of the 
Jacobian (2.19) can be solved efficiently. 

The technique of multiple shooting gives a robust method for determining long 
period UPOs in discrete dynamical systems. The idea also proves useful for continu- 
ous systems. Here there are several approaches one may take, for example, one may 
define a sequence of Poincare surface of sections (PSS) in such away that an orbit 
leaving one section reaches the next one in a predictable manner, without traversing 
other sections along the way, thus reducing the fiow to a set of maps. However, 
the topology of high- dimensional fiows is hard to visualise and such a sequence of 
PSSs will usually be hard to construct; see [17] for further details. We present an 
alternative approach when we discuss the application of Newton's method to fiows 
in 52.2.3. 
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The damped Newton-Raphson method 

The global convergence properties of the NR method are wildly unpredictable, and in 
general the basins of attraction cannot be expected to be simply connected regions. 
To improve upon this situation, it is useful to introduce the following function 

h{x) = ]^g{x)-g{x), (2.21) 

proportional to the square of the norm of g{x). Noting that the NR correction, 5x, 
is a direction of descent for /i, i.e. 

Vh-5x = {g''Dg)-{-Dg-'g), (2.22) 
= -/■^7<0, (2.23) 

suggests the following damped NR scheme 

Xi+i = Xi- Xi[Dg{xi)]-^g{xi), (2.24) 

where < Aj < 1 are chosen such that at each step h{x) decreases. The existence of 
such a Aj follows from the fact that the NR correction is in the direction of decreasing 
norm. 

At this stage it is important to make the distinction between the Newton step 
5x and the Newton direction d = —[Dg{x)]~^g{x). In the following, the Newton 
step will be a positive multiple of the Newton direction. The full Newton step 
corresponds to Aj = 1 in the above equation. In practice one starts by taking the 
full NR step since this will lead to quadratic convergence close to the root. In the 
case that the norm increases we backtrack along the Newton direction until we find 
a value of Aj for which h{xi+i) < h{xi). 

One such implementation, is given by the Newton-Armijo rule, which expresses 
the damping factor, Aj, in the form 

A, = 2-^ A; = 0,1,..., 

during each iteration successive values of k are tested until a value of k is found such 
that 

h{xi + \,d) < (1 - a\ifh{xi), (2.25) 
where the parameter a G (0, 1) is a small number intended to make (2.25) as easy 
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as possible to satisfy; see [53] for details and references concerning the Armijo rule. 

Methods like the Armijo rule are often called line searches due to the fact that 
one searches for a decrease in the norm along the line segment [xi ,Xi + d]. In 
practice some problems can respond well to one or two decreases in the step length 
by a modest amount (such as 1/2), whilst others require a more drastic reduction 
in the step-length. To address such issues, more sophisticated line searches can be 
constructed. The strategy for a practical line search routine is as follows: if after 
two reductions by halving, the norm still does not decrease sufficiently, we define 
the function 

$(A) := + Arf). (2.26) 

We may build a quadratic polynomial model of $ based upon interpolation at the 
three most recent values of A. The next A is the minimiser of the quadratic model. 
For example, let Aq = 1, Ai = 1/2 and A2 = 1/4, then we can model $(A) by 

<I>(A) ^ (^<l>o-8<l>i + y<l>2)A2 + (-2<l>o + 10$i-8<l>2)A + (^<l>o-2<l>i + ^<l>2), (2.27) 

where $j = $(Aj) for i = 0, 1,2. Taking the derivative of this quadratic, we find 
that it has a minimum when 

6<l>o - 30$! + 24$2 

16<l>o -48$! + 32$2' ' 

These ideas can of course be further generalised through the use of higher-order 
polynomials in the modeling of the function $(A). 

Damped NR is a very powerful detection routine. By using information contained 
in the norm of g to determine the correct step-size, a globally convergent method 
is obtained. Further, close to the root it reverts to the full NR step thus recovering 
quadratic convergence. It should, however, be pointed out that, in practice, the 
method may occasionally get stuck in a local minimum, in which case, detection 
should be restarted from a new seed. 



Quasi-Newton 

For an n-dimensional system, calculation of the Jacobian matrix requires partial 
derivative evaluations whilst approximating the Jacobian matrix by finite differences 
requires O(n^) function evaluations. Often, a more computationally efficient method 
is desired. 

Consider the case of the one-dimensional NR algorithm for solving the scalar 
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equation g{x) = 0. Recall that the secant method [77] approximates the derivative 
g'{xi) with the difference quotient 

k = (2.29) 

and then takes the step 

Xi+i = Xi- hr^g{xi). (2.30) 

The extension to higher dimensions is made by using the basic property of the 
Jacobian Dg{x)6x = 6g to write the equation 

Dg{xi){xi - Xi-i) = g{xi) - g{x^-i). (2.31) 

For scalar equations, (2.29) and (2.31) are equivalent. For equations in more than 
one variable, (2.29) is meaningless, so a wide variety of methods that satisfy the 
secant condition (or quasi-Newton condition) (2.31) have been designed. 

Broyden's method is the simplest example of the quasi-Newton methods. In the 
case of Broyden's method, if Xi and are the current approximate solution and 
Jacobian, respectively, then 

Xi+i = Xi- \iB~^g{xi), (2.32) 

where Aj is the step length for the approximate Newton direction di = —B^^g{xi). 
After each iteration, Bi is updated to form -Bi+i using the Broyden update 

(v — Bjs)s'^ 
s' s 

Here y = g{xi-^-i) — g{xi) and s = Xidi. It is a straightforward calculation to show 
that the above update formula satisfies the quasi-Newton condition (2.31). Different 
choices for the update formula are possible, indeed, depending on the form of the 
update a different quasi-Newton scheme results; see [53] for further details and 
references. 

2.2.3 Newton's method for flows 

Given an autonomous flow 

^=v{x), (2.34) 
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we can write the periodic orbit condition as 

^'^{x)-x = 0. (2.35) 

Here the flow map = x{t) corresponds to the solution of (2.34) at time t, and 
T defines the period. In order to determine UPOs of the flow we must determine 
the {n + 1) vector (x, T) satisfying Eq. (2.35). Immediately however, we run into 
a problem; namely, that the system in (2.35) is under determined. To solve this 
problem we add an equation of constraint by way of a Poincare surface of section 
(PSS). As long as the flow crosses the PSS transversely everywhere in phase space, 
we obtain a (n — l)-dimensional system of equations of full rank. In what follows 
we assume for simplicity that the PSS is given by a hyperplane V, that is, it takes 
the following linear form 

(.{x-xo) = 0, (2.36) 

where xq € P and is a vector perpendicular to V. The action of the constraint 
given by Eq. (2.36), is to reduce the n-dimensional flow to a (ri — l)-dimensional map 
defined by successive points of directed intersection of the flow with the hyperplane 
V. 

We may now proceed in similar fashion to the discrete case by linearising Eq. (2.35), 
we then obtain an improved guess by solving the resulting linear system in conjunc- 
tion with Eq. (2.36). The NR correction [6x,6T) is given by the solution of the 
following system of linear equations 



J -In v{<f>^{x)) 




Sx 




(f)^{x) — X 


C 




6T 








(2.37) 



where the corresponding Jacobian, J, is obtained by simultaneous integration of the 
fiow and the variational form [17] 

^ = AJ, J(0) = In, (2.38) 

here A = dv/dx and as usual I„ denotes the n x n identity matrix. Note that 
the second term in the vector on the RHS of Eq. (2.37) has been equated to zero, 
this is since our initial guess for the NR method lies on the PSS. By construction, 
Eq. (2.36) is equal to zero on the PSS. 

Since the ODE (2.34) is autonomous, it follows that if x{t) is a solution of (2.34), 
(2.35), then any phase shifted function x{t + g), g G M, also solves (2.34), (2.35). It 
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Figure 2.3: After integration time the point Xi returns to the Poincare surface of 
section (PSS), however the nearby point Xi + 5x does not. Thus, the matrix needed 
to map an arbitrary deviation 5xi on the PSS to the subsequent one fej+i needs to 
take into account the imphcit dependance of the return time on x. 



is this arbitrariness which forces us to introduce a PSS, by restricting to the surface 
V we ehminate corrections along the flow. Constraints can of course arise in other 
more natural ways, for example, if the flow has an integral of motion, such as the 
energy in a Hamiltonian system. In this case one may define a map of dimension 
(n — 2), by successive points of directed intersection of the flow, with the (n — 2)- 
manifold given by the intersection of the plane, P, and the corresponding energy 
shell. This leads to the following form for the NR correction 



J-I„ v{(f{x)) VH 
^"^0 



5x 
5T 
6E 



[X 



X 



(2.39) 



which must be solved simultaneously with 



H(x)~E = 0, 



(2.40) 



in order to remain on the energy shell. Here H is the Hamiltonian function and E 
defines its energy. 

To complete we mention that the preceding discussion seems somewhat strange, 
in that, in order to detect fixed points of the Poincare map which has dimension {n — 
1), we solve a system of {n + 1) equations. This is because, in general, the equations 
of constraint are nonlinear, recall that we chose to work with linear equations to 
simplify the exposition. In the case when the equations of constraint can be solved 
explicitly, one may work directly with the Poincare map, however, care must be taken 
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when computing the corresponding Jacobians, in particular, the total derivative 
should be taken due to the dependence of the constrained variables and the return 
times on x; see Figure 2.3. This may naturally be extended to the Hamiltonian case 
given above, and in fact to a flow with an arbitrary number of constraints; see [25] 
for further details. 

Multiple shooting revisited 

As in the discrete case, the detection of long period UPOs can be aided by the 
introduction of multiple shooting. We have seen that the problem of detecting UPOs 
of a flow may be regarded as a two-point boundary value problem; the boundary 
conditions being the relations closing the trajectory in phase space after time T, i.e. 
the periodic orbit condition, = x. In order to apply multiple shooting, we 

begin by discretising the time evolution into time intervals 

= To < ri < ■ ■ ■ < r^_i < rjv = T. (2.41) 

For ease of exposition we partition the time interval [0,T] into N intervals of equal 
length, A = T/N. The aim of multiple shooting is then to integrate N trajectories 
- one for each interval - and to check if the final values coincide with the initial 
conditions of the next one up to some precision. If not, then we apply NR to update 
the initial conditions. 

Let us denote by Xi the initial value of the trajectory at time Xj, and the final 
value of the trajectory at time Tj+i by 0'^(xj), then (A^ — 1) continuity conditions 
exist 

Ci{xi,Xi+i,T) = ^'^{xi)-Xi+, = 0, z = 0,1,..., AT -2, (2.42) 
together with the boundary conditions 

B{xn-i,Xo, T) = (p^ixN-i) -Xo = 0. (2.43) 

Now, we must solve A^ initial value problems, and for that we adopt the NR 
method 

a{x„ x.+i, T) + ^6x, + ^^S^i+i + = 0. (2.44) 

These equations become 



Ci{xi, Xi+i, T) + Ji5xi + 5xi+i + —Vi+i5T = 0, 



(2.45) 
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and using the boundary conditions (2.43), we get 



B{xn-i,Xo,T) + Jn-i5xn-i - 5xq + j^vnST = 0. 



Here 



J, 



dxi 



and Vi+i = v{(j)^{xi)). 



(2.46) 



(2.47) 



We can then write Eqs. (2.45), (2.46) in a matrix form of dimension [nN x n{N + 1)] 



Jo 

Ji -I, 



-In 













^^1 

V2 



Jn~2 —In VN-1 
Jn-1 Vn 



6xo 
Sxi 

6xn- 



N-2 



6xn-i 
—5T 



Co 
B 



which must be solved simultaneously with the PSS condition 



(2.48) 



P{x]Xo) 



0. 



(2.49) 



As in the case of simple shooting, the equation (2.49) ensures the Newton correction 
lies on the PSS. 

For a more detailed description in applying multiple shooting algorithms to flows 
- in particular conservative ones - see [26]. We conclude by mentioning the POMULT 
software described in [26] which is written in Fortran and provides routines for a 
rather general analysis of dynamical systems. The package was designed specifically 
for locating UPOs and steady states of Hamiltonian systems by using two-point 
boundary value solvers which are based on similar ideas to the multiple shooting 
algorithms discussed here. However, it also includes tools for calculating power 
spectra with fast fourier transform, maximum Lyapunov exponents, the classical 
correlation function, and the construction of PSS. 



2.3 Least-square optimisation tools 

In §2.2.2 we saw that by using the information contained within the norm of g{x) = 
f^{x) —X it is possible to greatly improve the convergence properties of the Newton- 
Raphson (NR) method. Let us now rewrite the function h{x) defined in Eq. (2.21) 
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in a slightly different form 

1 " 

M^) = 2E[^^(^)]'- (2-50) 

Here the gj are the components of the vector valued function g = f"^{x) — x. Since h 
is a nonnegative function for all x, it follows that every minimum x* of h satisfying 
h{x*) = is a zero of g, i.e. g{x*) = /^(x*) — a;* = 0. 

In order to detect UPOs of g we can apply the NR algorithm to the gradient of 
(2.50), the corresponding Newton direction is given by 

= -H-\xi)Vh{xi), (2.51) 

where H is the Hessian matrix of mixed partial derivatives. Note that the com- 
ponents of the Hessian matrix depend on both the first derivatives and the second 
derivatives of the gf 



Hki = J2 



99j dgj , . dgj 

OXi. OX] OXkOXi 



(2.52) 



In practice the second order derivatives can be ignored. To motivate this, note that 
the term multiplying the second derivative in Eq. (2.52) is gj{x) = fj{x) —Xj. For a 
sufficiently good initial guess the gj will be small enough that the term involving only 
second order derivatives becomes negligible. Dropping the second order derivatives 
leads to the Gauss-Newton direction 

di = -H-\xi)Vh{xi). (2.53) 

Here H = Dg'^Dg where Dg is the Jacobian of g evaluated at Xj. 

In the case where a good initial guess is not available, the application of Gauss- 
Newton (GN) will in general be unsuccessful. One way to remedy this is to move in 
the direction given by the gradient whenever the GN correction acts as to increase 
the norm, i.e. use steepest descent far from the root. 

Based on this observation, Levenberg proposed the following algorithm [59] . The 
update rule is a blend of the aforementioned algorithms and is given by 



Xi 



[H + Xln]-'Vh{xi). (2.54) 



After each step the error is checked, if it goes down, then A is decreased, increasing 
the influence of the GN direction. Otherwise, A is increased and the gradient dom- 
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inates. The above algorithm has the disadvantage that for increasingly large A the 
matrix (if + AI„) becomes more and more diagonally dominant, thus the information 
contained within the calculated Hessian is effectively ignored. To remedy this, Mar- 
quardt had the idea that an advantage could still be gained by using the curvature 
information contained within the Hessian matrix in order to scale the components 
of the gradient. With this in mind he proposed the following update formula 

Xi+i = x,-[H + X- diag{H)]-^Vh{xi), (2.55) 

known as the Levenberg-Marquardt (LM) algorithm [63]. Since the Hessian is pro- 
portional to the curvature of the function h, Eq. (2.55) implies larger steps in the 
direction of low curvature, and smaller steps in the direction with high curvature. 

LM is most commonly known for its application to nonlinear least square prob- 
lems for modeling data sets. However, it has recently been applied in the context of 
UPO detection. Lopez et al [60] used it to detect both periodic and relative periodic 
orbits of the complex Ginzburg-Landau equation. In this work Lopez et al used 
the Imder implementation of the LM algorithm contained in the MINPACK pack- 
age [66] . They found that in the case of the complex Ginzburg-Landau equation the 
method out-performs the quasi-Newton methods discussed earlier, as well as other 
routines in the MINPACK package such as hybr j which is based on a modification 
of Powell's hybrid method [76]. 

2.4 Variational methods 

For those equations whose dynamics are governed by a variational principle, periodic 
orbits are naturally detected by locating extrema of a so called cost function. In 
the case of a classical mechanical system such a cost function is given by the action, 
which is the time integral of the Lagrangian, L(g, g, t), where q G M". Starting from 
an initial loop, that is, a smooth closed curve, g(t), satisfying q{t) — q{t -|- T) = 0, 
one proceeds to detect UPOs by searching for extrema of the action functional 

S{q]= I Liq,q,t)dt. (2.56) 

One advantage of this approach is that UPOs with a certain topology can be detected 
since the initial guess is not a single point but a whole orbit. This method has 
recently been applied to the classical ri-body problem, where new families of solutions 



2.4 Variational methods 



36 



have been detected [89]. 

In the case of a general flow - such as that defined by Eq. (2.34) - one would 
still like to take advantage of the robustness offered by replacing an initial point by 
a rough guess of the whole UPO. To this end, Civtanovic et al [57] have introduced 
the Newton descent method, which can be viewed as a generalisation of the multiple 
shooting algorithms discussed earlier. The basic idea is to make an informed guess 
of what the desired UPO looks like globally, and then to use a variational method 
in order to drive the initial loop, L, towards a true UPO. 

To begin, one selects an initial loop L, a smooth, differentiable closed curve x{s) 
in phase space, where s G [0, 27r] is some loop parameter. Assuming L to be close to 
a UPO, one may pick (A^ — 1) pairs of nearby points along the loop and the orbit 

Xi = x{si), < Si < • ■ ■ < SAf_i < 2-71, (2-57) 
Xi = x{ti), < ti < ■■• < tjv-i < T. (2.58) 

Denote by 5xi the deviation of a point Xj on the periodic orbit from the point Xi. 
Let us define the loop velocity field, as the set of s-velocity vectors tangent to the 
loop L, i.e. v{x) = dx/ds. Then the goal is to continuously deform the loop L 
until the directions of the loop velocity field are aligned with the vector field v 
evaluated on the UPO. Note that the magnitudes of the tangent vectors depend 
upon the particular choice of parameter s. Thus, to match the magnitudes, a local 
time scaling is introduced 

\{si) = Ati/Asi. (2.59) 
Now, since Xi = Xi + 6xi lies on the periodic orbit, we have 

0At.+<5t^~. ^ ^ ^ (2.60) 

Linearisation of (2.60) yields the multiple shooting NR method 

Sxi+i - J{xi, Ati)Sxi - Vi+iSti = 0^*' (xi) - Xi+i, (2.61) 

which, for a sufficiently good guess, generates a sequence of loops L with a decreasing 
cost function given by 

F\x) = —— 5^(0'^*'(i.) - i.+i)^. (2.62) 
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As with any NR method, a decrease in the cost function is not always guaranteed if 
the full NR-step is taken. However, a decrease in F'^ is ensured if infinitesimal steps 
are taken. 

Fixing Asj one may proceed by (5r, where r is a fictitious time which parame- 
terises the Newton Descent, thus the RHS of Eq. (2.61) is multiphed by 5t 

5xi+i - J{xi, Ati)Sxi - Vi+i6ti = 5T{(j)^^'{xi) - Xi+i). (2.63) 

According to Eq. (2.59) the corresponding change in Ati is given by 

5ti = Ati{T + 6T)-Ati{T), 

dAU , 
-dr. 



dr 
dX 

= Asi—6T, (2.64) 
or 

and similarly 6x = {d/dT)x{si, t)6t. Dividing Eq. (2.63) by 6r and substituting the 
above expressions for 6ti, 6xi yields 

- Jixi, Ati)^ - v,+,^{si, T)Asi = 0^*'(xi) - Xi+i. (2.65) 
dr dr or 

Now in the limit N ^ oo, the step sizes Asj, Ati = 0{j^) ^0, giving 

Vi+i ~ Vi, Xj+i ~ Xj + ViAsi, 

J{xi, Ati) ~ / + A{x,)Ati, (p^'^ [xi] ^Xi + ViAti. 

Substituting into Eq. (2.65) and using the relation (2.59) results in the following 
PDE which describes the evolution of a loop L{t) toward a UPO 

d'^x , ,dx dX , ^ 

XA— - V— = Xv-v, (2.66) 



dsdr dr dr 

where v is the vector field given by the flow, and "v is the loop velocity field. The 
importance of Eq. (2.66) becomes transparent if we rewrite the equation in the 
following form 

d 

- Xv) = -[v - Xv)^ (2.67) 

or 
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from which it follows that the fictitious time flow decreases the cost functional 



The Newton descent method is an infinitesimal step variant of the damped NR 
method. It uses a large number of initial points in phase space as a seed which allows 
for UPOs of a certain topology to be detected, in practice it is very robust and has 
been applied in particular to the Kuramoto-Sivashinsky equation in the weakly 
turbulent regime, where many UPOs have been detected. The main problem - as 
with most multiple shooting methods - is that of efficiency. In [57] they state that 
most of the computational effort goes into inverting the [{nd + 1) x {nd + 1)] matrix 
resulting from the numerical approximation of Eq. (2.66), here d is the number of 
loop points and will typically be large. 

2.5 Summary 

Our review is not exhaustive and further methods exist for the detection of unstable 
periodic orbits (UPO) in a chaotic dynamical system. For example, in those cases 
where the systems dynamics depend smoothly upon some parameter the method of 
continuation can be applied. Here one uses the fact that their exists a window in 
parameter space where one can easily detect solutions, the idea is then to slowly 
vary the parameter into a previously unaccessible region in parameter space whilst 
carefully following the corresponding solution branch. These ideas have recently 
been applied to determine relative periodic motions to the classical fluid mechanics 
problem of pressure-driven flow through a circular pipe [97]. Another recent proposal 
by Parsopoulos and Vrahatis [71] uses the method of particle swarm optimization 
in order to detect UPOs as global minima of a properly defined objective function. 
The preceding algorithm has been tested on several nonlinear mappings including 
a system of two coupled Henon maps and the Predator-Prey mapping, and has 
been found to be an efficient alternative for computing UPOs in low-dimensional 
nonlinear mappings. 

For generic maps (or flows) no computational algorithm is guaranteed to flnd 
all solutions up to some period p (time T^ax)- For systems where the topology is 
well understood one can build sufficiently good initial guesses so that the Newton- 
Raphson method can be applied successfully. However, in high-dimensional prob- 
lems the topology is hard to visualize and typically Newton-Raphson will fail. Here 




(2.68) 
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the variational method offers a robust alternative, methods that start with a large 
number of initial guesses along an orbit, do not suffer from the numerical instability 
caused by exponential sensitivity of chaotic trajectories. Unfortunately, these meth- 
ods tend to result in large systems of equations, the solution of which is a costly 
process, this has the knock on effect of slowing the convergence of such methods 
considerably. This makes it extremely difficult to find UPOs with larger periods or 
detect them for higher dimensional systems. In Chapters 3 and 4 we shall introduce 
and discuss the method of stabilising transformations in some detail, due to the fact 
that the method possesses excellent global convergence properties and needs only 
marginal a priori knowledge of the system, it is ideally suited to determining UPOs 
in high- dimensional systems. 



Chapter III 



Stabilising transformations 



Tis plain that there is not in nature a point of stabihty to be found: 
everything either ascends or dechnes. 
Sir Walter Scott 

An algorithm for detecting unstable periodic orbits based on stabilising transfor- 
mations (ST) has had considerable success in low- dimensional chaotic systems [19]. 
Applying the same ideas in higher dimensions is not trivial due to a rapidly in- 
creasing number of required transformations. In this chapter we introduce the idea 
behind the method of STs before going onto to analyse their properties. We then 
propose an alternative approach for constructing a smaller set of transformations. 
The performance of the new approach is illustrated on the four-dimensional kicked 
double rotor map and the six-dimensional system of three coupled Henon maps [14]. 

3.1 Stabilising transformations as a tool for de- 
tecting UPOs 

The method of STs was first introduced in 1997 by Schmelcher and Diakonos (SD) as 
a tool for determining unstable periodic orbits (UPOs) for general chaotic maps [83]. 
For the first time UPOs of high periods where determined for complicated two- 
dimensional maps, for example, the Henon map [43] and the Ikeda-Hammel-Jones- 
Maloney map [50, 38]. One of the main advantages of the method is that the basin 
of attraction for each UPO extends far beyond its linear neighbourhood and is given 
by a simply connected region. 
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Let us consider the following n-dimensional system: 



U: Xi+i 



fix,), /:R"^M^ 



(3.1) 



The basic idea of the SD method is to introduce a new dynamical system U, with 
fixed points in exactly the same position in phase space as the UPOs of U but with 
differing stability properties, ideally the fixed points of U will be asymptotically 
stable. In general however, no such system exists. Fortunately, this can be rectified 
by choosing a set of associated systems, such that, for each UPO of U, there is at 
least one system which possesses the UPO as a stable fixed point. To this end, SD 
have put forward the following set of systems 



here A is a small positive number, p is the period and Ck is an n x n matrix with 
elements [Ck]ij G {0, ±1} such that each row or column contains only one nonzero 
element. To motivate Eq. (3.2), it is useful to look at the Jacobian of the system Uk 



For a fixed point of Uk to be stable all eigenvalues of the above expression must have 
absolute value less than one. This can be achieved in the following way: (i) choose 
a linear transformation Ck such that all eigenvalues of the matrix Ck[DfP{xi) — I„] 
have negative real part, and (ii) pick A sufficiently small so that the eigenvalues of 
I + XCklDf^lxi) — In] are scaled so as to all have absolute value less than one. The 
fact that such a A exists is clear, whilst step (i) follows from Conjecture 1.1. 

This leads to the following algorithm to detect all period-p orbits of a chaotic 
map. We begin by placing seeds over the attractor using any standard seeding 
strategy, for example, a chaotic trajectory or a uniform grid. Using the different 
matrices Ck, we construct the 2"n! associated systems Uk - see Observation 1.1. For 
a sufficiently small value of A, we propagate each seed using the iteration scheme 
in (3.2) to compute a sequence {xi} for each of the 2"n! maps, Uk- If a sequence 
converges, we check whether a new period-p orbit point has been found, and if so, 
proceed to detect the remaining [p — 1) orbit points by iterating the map /. In 
order to test for the completeness of a set of period-p orbits, SD suggest iterating a 
multiple of the initial seeds used in the detection process. If no new orbits are found 
then they claim that the completion of the set is assured. 



Uk- Xi+i = Xi + XCkirixi) - Xi], k = 1,. . . , 2"n!, 



(3.2) 



(ixj_|_i 



ln + \Ck[Dr{x,)-ln]. 



(3.3) 
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From a computational perspective the SD method has two major failings. Firstly, 
and most importantly, the application to high- dimensional systems is restricted due 
to the fact that the number of matrices in the set Csd increases very rapidly with 
system dimension. Secondly, it follows from our prior discussion that the parameter 
A scales with the magnitude of the largest eigenvalue of (3.3). Since the instability 
of a UPO increases with the period p, very small values of A must be taken in order 
to detect long period orbits. This has the effect of increasing the number of steps - 
hence function evaluations - needed to obtain convergence, which in turn leads to 
a rise in the computational costs. 

The problem of extending the method of STs to high-dimensional systems is the 
primary concern of this thesis and will be dealt with in some detail in this chapter 
and the next. The problem of efficiency is solved to some extent by the following 
observation in [84]. 

Observation 3.1. For a given map f. Taking the limit X ^ in Eq. (3.2) leads 
to the following continuous flow 

lim ^^±i-^ = ^ = anx) -x) = C,g{x). (3.4) 
\~*o A as 

The equilibria of Eq. (3.4) are located at exactly the same positions and share the 
same stability properties as the fixed points of the associated systems Uk in the 
limit of small A. By transferring to the continuous setting, the dependency on A has 
been removed, this allows us to use our favourite off-the-shelf numerical integrator 
to enable us to detect UPOs of /; the SD method is precisely Euler's method to 
solve the fiow of (3.4) with step-size At = A. 

Taking into account the typical stiffness of the flow in Eq. (3.4), Davidchack 
and Lai proposed a modifled scheme employing a semi-implicit Euler method [19]. 
The semi-implicit scheme is obtained from the fully implicit Euler method in the 
following way: starting from the implicit scheme as applied to Eq. (3.4) 

Xj+i = Xi + hCg{xi+i), (3.5) 

one obtains the semi- implicit Euler routine by expanding the term (^(xj+i) in a 
Taylor series about Xi and retaining only those terms which are linear 



Xj+i = Xi + hCg{xi) + hCg\xi)Axi + 0{Ax1)., 
= Xi + [\c'^ - g\xi)Y^g{xi). 



(3.6) 
(3.7) 
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Figure 3.1: (colour online) Shown in red are the basins of convergence of (a) the 
Newton method, (b) the Schmelcher and Diakonos method with < A < 0.3568 and 
C = 1, and (c) the Davidchack and Lai method with = 4.0 and C = I to the zeros 
of a function g{x) = cos(x^) in the interval (—3, 3). Arrows indicate the direction of 
convergence, and large dots are the zeros to which the method converges. 

The Davidchack-Lai method is completed by choosing the step-size h = 

here /? > is a scalar parameter and Sj = [[(^(xj)!! is an L2 norm. This leads to the 

following algorithm for detecting UPOs 

=Xi + [Ps.C - G,]-'gix,), (3.8) 

where Gi = Dg{xi) is the Jacobian matrix, and "T" denotes transpose. 

In the vicinity of an UPO, the function g{x) tends to zero and Newton-Raphson is 
restored, indeed it can be shown that the method retains quadratic convergence [55]. 
Away from the UPO and for sufficiently large (3, the flow is accurately reproduced 
thus conserving the global convergence properties of the SD method. A nice property 
of the above algorithm is that the basin sizes can be controlled via the parameter (3. 
In particular in [21] it is shown that increasing the value of the parameter (3 results 
in a larger basin size. A large value of (3, however, requires more time steps for 
the iteration to converge to a UPO. This leads to a trade off between enlarging the 
basins and the speed of convergence. A comparison of the basins of attraction for 
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the Newton-Raphson, the Schmelcher-Diakonos and the Davidchack-Lai methods is 
given in Figure 3.1. Here the roots of the function g{x) = cos(x^) are shown on 
the interval (—3, 3) along with the corresponding basins of attraction for the three 
methods; see [19] for further details. 

3.1.1 Seeding with periodic orbits 

An important ingredient of the Davidchack-Lai algorithm lies in the selection of 
initial seeds. Davidchack and Lai claim that the most efficient strategy for detecting 
UPOs of period p is to use UPOs of other periods as seeds [19]. This can be 
understood due to the distribution of the UPOs. It is well known that orbit points 
cover the attractor in a systematic manner, which in turn reflects the foliation of 
the function f^{x) and its iterates. For low-dimensional maps, Davidchack and Lai 
have considerable success with the aforementioned seeding strategy as compared to 
traditional seeding algorithms. Indeed, for the Henon and the Ikeda-Hammel-Jones- 
Maloney maps, a//period-p orbits were detected using period-(p — 1) orbits, provided 
they exist; in the case that no period- (p — 1) orbits exist, one can try either {p — 2) 
or {p + 1). 

The use of UPOs as seeds is also crucial to our work. We shall see that it is the 
information contained in the UPOs - or rather their Jacobian's - that will enable 
us to construct new sets of STs, in turn allowing for efficient detection of UPOs in 
high-dimensional systems. 

3.2 Stabilising transformations in two dimensions 

The stability of a fixed point x* of the flow 

E: ^ = Cgix) , (3.9) 

is determined by the real parts of the eigenvalues of the matrix CG, where G = 
Dg{x*) is the Jacobian matrix of g{x) evaluated at x*. For x* to be a stable flxed 
point of E, the matrix C has to be such that all the eigenvalues of CG have negative 
real parts. In order to understand what properties of G determine the choice of a 
particular stabilising transformation C, we use the following parametrisation for the 
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general two-dimensional orthogonal matrices 



s COS a sm a 
—s sin a cos a 



(3.10) 



where s = ±1 and — vr < a < vr. When a = —it/2, 0, tt/2, or vr, we obtain the set 
of matrices Csd- For example, 



and 



If we write 



Cl,7r/2 



C-l,7r 



1 
-1 

1 
-1 



G:={G,j}, I, J = 1,2, 
then the eigenvalues of Cs,aG are given by the following equations: 



(3.11) 



(3.12) 



(3.13) 



cri_2 = -Acos(a - 9) ± ^/A'^cos'^{a -9) - sdetC 



(3.14) 



where det G = G11G22 - ^12^21, A = | ^(sGn + ^22)^ + (sGi2 - G21)', and 

sGi2 — ^21 



tan^ 



— sGii — G 



-Tl < 9 < Tl . 



(3.15) 



22 



Note that the signs of the numerator and denominator are significant for defining 
angle 9 in the specified range and should not be canceled out. It is clear from 
Eq.(3.14) that both eigenvalues have negative real parts when 



s = s := sgn det G, and |a; — 6'| < f 



2 ' 



(3.16) 



provided that det G 7^ 0. This result proves the validity of Conjecture 1.1 for n = 2. 
Moreover, it shows that there are typically two matrices in Csd that stabilise a given 
fixed point. 

Parameter 9 clearly plays an important role in the above analysis. The following 
two theorems show its relationship to the eigenvalues and eigenvectors of the stability 
matrix of a fixed point. 

Theorem 3.1. Let x* be a saddle fixed point of f^{x) : 1-^ whose stability 
matrix Df^{x*) has eigenvalues Ai_2 such that IA2I < 1 < |Ai| and eigenvectors 
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defined by the polar angles < 0i_2 < ^r, i.e. Vi^2 = (cos 0i_25 sin 01^2)"'" ■ Then the 
following is true for the angle 9 defined in Eq. (3.15): 



Case 1. Ai < -1: 



Moreover, if ^ 1, then 

^ = (0i-02)(mod vr) - f + . 

Case 2. Ai > 1: 



e 



37r 
2 



/>i - </>2 , <</>!- </)2 < vr , 

>1 - 02 , -TT < 01 - 02 < . 



(3.17) 



(3.18) 



(3.19) 



Proof. Matrix G = Df^{x*) — /2, where I2 is the identity matrix, can be written as 
follows: 



G 



Gil 


G12 




G21 


G22 





cos 01 COS 02 

sin 01 sin 02 



Ai - 1 
A2- 1 



cos 01 COS 02 

sin 01 sin 02 



(3.20) 



Case 1: Since det G = (Ai — 1)(A2 — 1) > we set s = 1 and obtain from Eq. (3.15): 

(Ai - A2) COt(0i - 02) 



tan^ 



2 - Ai - A2 



(3.21) 



where, just like in Eq. (3.15), as well as in Eqs. (3.22) and (3.25) below, the 
signs of the numerator and denominator should not be canceled out. Since 
2 — Ai — A2 > 0, we have that cos 6' > or 



0G(-f,f) . 



For I All > 1, Eq. (3.21) yields: 



tan^ 



All 



cot (01 — 



1 - V2 



and, given Eq. (3.17), the result in Eq. (3.18) immediately follows. 
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Case 2: In this case det G = (Ai — 1)(A2 — 1) < 0, so, from Eq. (3.15) with s = — 1: 

(A2 - Ai) cos(0i + 02)/ sin(0i - 02) 



(A2 - Ai) sin(0i + 02)/ sin(0i 
_ - cos(0i + 02)/ sin(0i - 02) 
- sin(0i + 02)/ sin(0i - 02) ' 

since A2 — Ai < 0. The result in Eq. (3.19) follows. 

Theorem 3.2. Let x* be a spiral fixed point of p{x) : ^ 
matrix Df^{x*) has eigenvalues Ai^2 = A ± icu. Then 

G (-f,f) if A<1, 

e e (-7r,-f) U (f,7r) if A>1. 

Proof. The stability matrix can be decomposed as follows: 
DfP{x*) 



(3.22) 



□ 



whose stability 



(3.23) 



cosq 


i e'' 




" A 


UJ 




cos 


e"" 


sin 4 


) 




— UJ 


A 




sin 






(3.24) 



where G M. Given that G = Df^{x*) — I2, we have from Eq. (3.15): 

-UJ cosh?]/ sin0 



tan 6' 



1-A 



The result in Eq. (3.23) follows from the sign of the denominator. 



(3.25) 
□ 



The key message of the above theorems is that the ST matrix depends mostly on 
the directions of the eigenvectors and the signs of the unstable^ eigenvalues of Df^ 
(or their real parts), and only marginally on the actual magnitudes of the eigenvalues. 
This means that a transformation that stabilises a given fixed point x* of will also 
stabilise fixed points of all periods with similar directions of eigenvectors and signs 
of the unstable eigenvalues. In the next Section, we will show how this observation 
can be used to construct STs for efficient detection of periodic orbits in systems 
with n > 2. 



"'^That is, eigenvalues whose magnitude is larger than one. 
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3.3 Extension to higher-dimensional systems 

To extend the analysis of the preceding Section to higher-dimensional systems, we 
note that the matrix Cgfi, as defined by Eqs. (3.10), (3.15), and (3.16), is closely 
related to the orthogonal part of the polar decomposition of G; see Appendix D. 
Recall that any non-singular nxn matrix can be uniquely represented as a product 



G = QB, 



(3.26) 



where Q is an orthogonal matrix and i? is a symmetric positive definite matrix. The 
following theorem provides the link between Cs^e and Q for n = 2: 

Theorem 3.3. Let G G M^^^ be a non-singular matrix with the polar decomposition 
G = QB, where Q is an orthogonal matrix and B is a symmetric positive definite 
matrix. Then matrix Gg^, as defined by Eqs. (3.10), (3.15) and (3.16), is related to 
Q as follows: 

Gs,e = -0" (3.27) 

Proof. Since Ggfi is an orthogonal matrix by definition, it is sufficient to prove that 
GgfiG is symmetric negative definite. Then, by the uniqueness of the polar decom- 
position, it must be equal to —B. 

Denote by 6jj the element in the i-th row and j-th column of G^^G. We must 
show that bi2 = &2i- Using Eq. (3.15), we have that 



n2 



and similarly 



sGi2 COS 9 + sin 9 



22 ' 



sGi2 + G22 — ~ 



G 



21 



-sGn 
GuGi2 + G21G22 



sGn + G' 



22 



- G22 
cos 9 . 



(3.28) 



cos^ 



J21 



G21 COS 9 — sGu sin 9 

sGi2 — G21 

(jr21 — SUrii 



-sGu 
G11G12 + G21G22 



sGu + G 



22 



- G22 
cos 9 . 



(3.29) 



cos 6* 



hence the matrix Gg^G is symmetric. Since, by definition, 9 and s are chosen such 
that the eigenvalues of Gg^G are negative, the matrix Gg^G is negative definite. 
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Finally, by the uniqueness of the polar decomposition, 

CsfiG = ~B = —Q^G , 

which completes the proof. □ 

For > 2, we can always use the polar decomposition to construct a trans- 
formation that will stabilise a given fixed point. Indeed, if a fixed point x* of an 
77,-dimensional flow has a non-singular matrix G = Dg{x*), then we can calculate 
the polar decomposition G = QB and use 

C = -Q^ , (3.30) 

to stabilise x* . Moreover, by analogy with the two-dimensional case, we can expect 
that the same matrix G will also stabilise fixed points x with the matrix G = Dg{x), 
as long as the orthogonal part Q of the polar decomposition G = QB is sufficiently 
close to Q. More precisely. 

Observation 3.2. G will stabilise x, if all eigenvalues of the product QQ^ have 
positive real parts. 

We base this observation on the following corollary of Lyapunov's stability the- 
orem; see Appendix D 

CoROLLORY 3.1. Let B G W"^"- he a positive definite symmetric matrix. If Q & 
R"'^"' is an orthogonal matrix such that all its eigenvalues have positive real parts, 
then all the eigenvalues of the product QB have positive real parts as well. 

Proof. According to Lyapunov's theorem, a matrix A G R"^" has all eigenvalues 
with positive real parts if and only if there exists a symmetric positive definite 
G e M"''" such that GA + A^G = H is positive definite. 

Let A = QB and let's choose G in the form G = \QB^^Q'^ . Since B is positive 
definite, its inverse B~^ is also positive definite, and, since G and B~^ are related 
by a congruence transformation, according to Sylvester's inertia law - see Appendix 
D - G is also positive definite. Now, 

GQB + {QBYG = IQB-'Q'^QB + ^BQ'^QB-'Q^ = i[g + Q^] . 

Therefore, QB has eigenvalues with positive real parts if and only if | [Q + Q'^] is 
positive definite. The proof is completed by observing that, for orthogonal matrices, 
the eigenvalues of + Q"""] are equal to the real parts of the eigenvalues of Q. □ 
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Note that Observation 3.2 is a direct generalisation of conditions in Eq. (3.16) 
which are equivalent to requiring that the eigenvalues of Cg^aCj^e have positive real 
parts. 

In the scheme where already detected periodic orbits are used as seeds to detect 
other orbits [19], we can use C in Eq. (3.30) as a stabilising matrix for the seed 
X*. Based on the analysis in §3.2, this will allow us to locate a periodic orbit in 
the neighbourhood of x* with similar invariant directions and the same signs of 
the unstable eigenvalues. Note, however, that the neighbourhood of the seed x* 
can also contain periodic orbits with the similar invariant directions but with some 
eigenvalues having the opposite sign (i.e. orbits with and without reflections). To 
construct transformations that would stabilise such periodic orbits, we can determine 
the eigenvalues and eigenvectors of the stability matrix of x* 

DfP{x*) = VAV-^ , (3.31) 

where A := diag(Ai, . . . , A„) is the diagonal matrix of eigenvalues of Df^{x*) and 
V is the matrix of eigenvectors, and then calculate the polar decomposition of the 
matrix 

G = V{SA-ln)V-\ (3.32) 

where 5* = diag(±l, ±1, . . . , ±1). Note that, as follows from the analysis in §3.2 for 
n = 2 and numerical evidence for n > 2, changing the sign of a stable eigenvalue 
will not result in a substantially different ST. Therefore, we restrict our attention 
to the following subset of S: 

Sii = I ' ^ *^ ' for i = l,...,n. (3.33) 
I 1, |A^|<1, 

For a seed with k real unstable eigenvalues, this results in 2^ possible transforma- 
tions. Note that, on the one hand, this set is much smaller than Csd, while, on 
the other hand, it allows us to target all possible types of periodic orbits that have 
invariant directions similar to those at the seed. 

3.4 Numerical results 

In this Section we illustrate the performance of the new STs on a four-dimensional 
kicked double rotor map [80] and a six- dimensional system of three coupled Henon 
maps [75]. Both systems are highly chaotic and the number of UPOs is expected to 
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grow rapidly with increasing period. The goal is to locate all UPOs of increasingly 
larger period. Of course, the completeness of the set of orbits for each period cannot 
be guaranteed, but it can be established with high degree of certainty by using the 
plausibility criteria outlined in the Introduction. 

In order to start the detection process, we need to have a small set of periodic 
orbits (of period p > 1) that can be used as seeds. Such orbits can be located 
using, for example, random seeds and the standard Newton-Raphson method (or 
the scheme in Eq. (3.8) with (3 = 0). We can then use these periodic orbits as seeds 
to construct the STs and detect more UPOs with higher efficiency. The process 
can be iterated until we find no more orbits of a given period. In previous work 
DL [19, 21] showed that for two-dimensional maps such as Henon and Ikeda it is 
sufficient to use period-(p— 1) orbits as seeds to locate plausibly all period-p orbits. 
For higher-dimensional systems, such as those considered in the present work, these 
seeds may not be sufficient. However, it is always possible to use more seeds by, for 
example, locating some of the period-(p+ 1) orbits, which can then be used as seeds 
to complete the detection of period-p orbits. The following recipe can be used as a 
general guideline for developing a specific detection scheme for a given system: 

1. Find a set of orbit points of low period using random seeds and the iterative 
scheme in Eq. (3.8) with /3 = (i.e. the Newton-Raphson scheme). 

2. To locate period-p orbits, first use period-{p — 1) orbits as seeds. For each 
seed xo, construct 1^ STs C using Eqs. (3.31-3.33), where k is the number of 
unstable eigenvalues of Df^^^{xo) . 

3. Starting from xq and with a fixed value of j3 > use the iteration scheme in 
Eq. (3.8) to construct a sequence {xj} for each of the 2^ STs. If a sequence 
converges to a point, check whether it is a new period-p orbit point, and if so, 
proceed to find a complete orbit by iterating the map f . 

4. Repeat steps 2 — 3 for several (3 in order to determine the optimal value of this 
parameter (see explanation below). 

5. Repeat steps 2 — 4 using newly found period-p points as seeds to search for 
period-{p + 1) orbits. 

6. Repeat steps 2 — 4 using incomplete set of period-{p + 1) orbits as seeds to find 
any missing period-p orbits. 
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Although we know that the action of (3 is to increase the basin size of the sta- 
bilised points, it is not known a priori what values of (3 to use for a given system 
and period. Monitoring the fraction of seeds that converge to periodic orbits, we ob- 
serve that it grows with increasing j3 until it reaches saturation, indicating that the 
iterative scheme faithfully follows the flow S. On the other hand, larger (3 translates 
into smaller integration steps and, therefore, longer iteration sequences. Thus the 
optimal value of [3 is just before the saturation point. As demonstrated previously 
by DL [19] and observed in the numerical examples presented in the following sec- 
tions, this value appears to scale exponentially with the period and can be estimated 
based on the information about the detection pattern at lower periods. 

The stopping criteria in step 3, which we use in the numerical examples discussed 
below, are as follows. The search for UPOs is conducted within a rectangular region 
containing a chaotic invariant set. The sequence {xj} is terminated if (i) Xj leaves 
the region, (ii) i becomes larger than a pre-defined maximum number of iterations 
(we use i > 100 + 5/3 ), (iii) the sequence converges, such that ||5f(a;j)|| < Tolg. In 
cases (i) and (ii) a new sequence is generated from a different seed and/or with a 
different stabilising matrix. In case (iii) five Newton iterations are applied to 
to allow convergence to a fixed point to within the round-off error. A point x* for 
which ||5f(x*)|| is the smallest is identified with a fixed point of The maximum 
round-off error over the set Xp of all detected period-p orbit points 

e^ax(p) = max{\\g{x*)\\ : x* E Xp} (3.34) 

is monitored in order to assess the accuracy of the detected orbits. 

To check if the newly detected orbit is different from those already detected, its 
distance to other orbit points is calculated: if ||x* — y*\\oQ > Tol^ for all previously 
detected orbit points y*, then x* is a new orbit point. Even for a large number 
of already detected UPOs, this check can be done very quickly by pre-sorting the 
detected orbit points along one of the system coordinates and performing a binary 
search for the points within Tol^ of x*. The infinity norm in the above expression 
is used for the computational efficiency of this check. 

The minimum distance between orbit points 

dminip) =min{||x* - y*\\oo ■ x*,y* E Xp] (3.35) 

is monitored and the algorithm is capable of locating all isolated UPOs of a given 
period p as long as emax(p) < Tolg < Tol^ < d^i^{p). Since typically emax(p) increases 
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and djainip) decreases with p (see Tables 3.1 and 3.2), the above conditions can be 
satisfied up to some period, after which higher-precision arithmetics needs to be used 
in the evaluation of the map. For the numerical examples presented in the following 
sections we use double-precision computation with Tolg = 10~^ and Tol^ = 10~^. 

3.4.1 Kicked double rotor map 

The kicked double rotor map describes the dynamics of a mechanical system known 
as the double rotor under the action of a periodic kick; a derivation is given in 
Appendix C. It is a 4-dimensional map defined by 

_yi+i 

where Xi G are the angle coordinates and yi G are the angular velocities after 
each kick. Parameters L and M are constant 2x2 matrices that depend on the 
masses, lengths of rotor arms, and friction at the pivots, while c G is a constant 
vector whose magnitude is proportional to the kicking strength /q. In our numerical 
tests we have used the same parameters as in [80], with the kicking strength /o = 8.0. 

The following example illustrates the stabilising properties of the transforma- 
tions constructed on the basis of periodic orbits. Let us take a typical period-3 orbit 
point X* = (0.6767947,5.8315697), y* = (0.9723920,-7.9998313) as a seed for lo- 
cating period-4 orbits. The Jacobian matrix Df^{x*,y*) of the seed has eigenvalues 
A = diag(206.48, -13.102, -0.000373,0.000122). Therefore, based on the scheme 
discussed in §3.3 Eqs. (3.31-3.33), we can construct four STs C corresponding to 
(5*11, 5*22) in Eq. (3.33) being equal to (+,+), (-,+), (+, -) and (-,-). Of the 
total of 2190 orbit points of period-4 (see Table 3.1), the transformations Ci, C2, 
C3, and C4 stabilise #(1) = 532, #(2) = 544, #(3) = 474, and #(4) = 516 orbit 
points, respectively, and these sets of orbits are almost completely non-overlapping. 
That is, the number of orbits stabilised by both Ci and C2 is #(ln2) = 2. Similarly, 
#(1 n 3) = 16, #(1 n 4) = 0, #(2 n 3) = 0, #(2 n 4) = 14, and #(3 n 4) = 0. On 
the other hand, the number of period-4 orbits stabilised by at least one of the four 
transformations is#(lU2U3U4)= 2034. This is a typical picture for other seeds 
of period-3 as well as other periods. 

This example provides evidence for the validity of our approach to constructing 
the STs in high-dimensional systems based on periodic orbits. It also shows that, 
in the case of the double rotor map, a single seed is sufficient for constructing 



Myi + Xi (mod 27r) 
Lyi + c sin Xi+i 



(3.36) 
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Table 3.1: Number n{p) of prime period-p UPOs, and the number N{p) of fixed 
points of p-times iterated map for the kicked double rotor map. The asterisk for 
p = 8 indicates that this set of orbits is not complete. Parameters eniax(p) and 
{p) are defined in Eqs. (3.34) and (3.35). 



p 


n{p) 


Nip) 


emax(p) 


dminip) 


1 


12 


12 


1.0- 


10-14 


1.3 


• 10" 


2 


45 


102 


5.9- 


10-14 


3.4- 


10-1 


3 


152 


468 


5.8- 


10-13 


6.2- 


10-2 


4 


522 


2190 


2.7- 


10-12 


6.9- 


10-3 


5 


2200 


11012 


2.6- 


10-11 


1.1 ■ 


10-3 


6 


9824 


59 502 


1.6- 


10-10 


1.8- 


10-^ 


7 


46 900 


328 312 


9.7- 


10-10 


9.1 ■ 


10-5 


8* 


229 082 


1834 566 


1.2 


•10-8 


5.5- 


10-5 



transformations that stabilise majority of the UPOs. Of course, in order to locate the 
UPOs, we need to ensure that the seeds are in the convergence basins of the stabilised 
periodic orbits. That is why we need to use more seeds to locate plausibly all periodic 
orbits of a given period. Still, because of the enlarged basins of the stabilised orbits, 
the number of seeds is much smaller than that required with iterative schemes that 
do not use the STs. 

Compared to the total of 384 matrices in Csd, we use only two or four transfor- 
mations for each seed, depending on the number of unstable directions of the seed 
orbit points. Yet, the application of the detection scheme outlined in §3.4 allows us 
to locate plausibly all periodic orbits of the double rotor map up to period 7. Table 
3.1 also includes the number of detected period-8 orbits that were used as seeds to 
complete the detection of period 7. 

The confidence with which we claim to have plausibly complete sets of periodic 
orbits for each period is enhanced by the symmetry consideration. That is, since the 
double rotor map is invariant under the change of variables {x,y) ^ (27r — x, —y), 
a necessary condition for the completeness of the set of orbits for each period is 
that for any orbit point {x*,y*) the set also contains an orbit point (2tt — x*, —y*)- 
Even though this condition was not used in the detection scheme, we find that the 
detected sets of orbits (apart from period 8) satisfy this symmetry condition. Of 
course, this condition is not sufficient to prove the completeness of the detected 
sets of UPOs, but, combined with the exhaustive search procedure presented above, 
provides a strong indication of the completeness. 
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3.4.2 Coupled Henon maps 

Another system we use to test the efficacy of our approach is a six-dimensional 
system of three coupled Henon maps (CHM) , 

xl^^ = a - {xiy + bxl_^, for j = 1,2,3, (3.37) 

where a = 1.4 and b = 0.3 are the standard parameter values of the Henon map and 
the coupling is given by 

xl = {l-e)xl + le{xl^'+xl-'), (3.38) 

with x^ = xf and xf = xj. We have chosen the coupling parameter e = 0.15. Our 
choice of this system is motivated by the work of Politi and Torcini [75] in which 
they locate periodic orbits in CHM for a small coupling parameter by extending 
the method of Biham and Wenzel (BW) [3]. This makes the CHM an excellent test 
system, since we can compare our results against those for the BW method. The 
BW method defines the following artificial dynamics 

xKt) = i-iy^^'^ixu,^ - a + - bxim, (3.39) 

with s{i,j) G {0, 1}. Given the boundary condition Xp_^i = x{, the equilibrium states 
of Eq. (3.39) are the period-p orbits for the CHM. The BW method is based on the 
property that every equilibrium state of Eq. (3.39) can be made stable by one of the 
possible sequences of s{i,j) and, therefore, can be located by simply integrating 
Eq. (3.39) to convergence starting from the same initial condition ocl = 0.0. It is 
also found that, for the vast majority of orbits, each orbit is stabilised by a unique 
sequence of s{i,j). 

In order to reduce the computational effort, Politi and Torcini suggest reducing 
the search to only those sequences s{i,j) which are allowed in the uncoupled system, 
i.e. with e = 0. This reduction is possible because the introduction of coupling has 
the effect of pruning some of the orbits found in the uncoupled Henon map without 
creating any new orbits. 

We have implemented the BW method with both the full search and the reduced 
search (BW-r) up to as high a period as is computationally feasible (see Table 3.2). 
In the case of the full search we detect UPOs up to period 8 and in the case of the 
reduced search up to period 12. The seed = 0.0 was used for all periods except 
for period 4, where it was found that with this seed both BW and BW-r located 
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Table 3.2: The number of prime UPOs for the system of three coupled Henon maps 
(CHM) detected by three different methods: BW - full Biham-Wenzel, BW-r - 
reduced Biham-Wenzel, ST - our method based on stabilising transformations, Max 
- maximum number of detected UPOs obtained from all three methods and the 
system symmetry. See text for details. 
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dmin(p) 


1 


8 


8 


8 


8 


1.3 • 




9.9 


10-1 


2 


28 


28 


28 
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10-3 
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64 
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10-^ 
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10-1 
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563 


568 
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10-1 
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7.1 


10-9 


5.4 


lo-'^ 


12 
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1999 
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2.5 


10-6 


4.3 


10-1 


13* 






1079 




8.6 


10-8 


4.0 


lo-'^ 


14* 






6599 




2.3 


10-6 


3.5 


lo-'' 


15* 






5899 




7.0 


10-6 


1.5 


lo-'' 



only 28 orbits. We found a maximum of 34 orbits using the seed xl = 0.5. It is 
possible that more orbits can be found with different seeds for other periods as well, 
but we have not investigated this. The example of period 4 illustrates that, unlike 
for a single Henon map, the Biham-Wenzel method fails to detect all orbits from a 
single seed. 

Even though our approach (labeled "ST" in Table 3.2) is general and does not 
rely on the special structure of the Henon map, its efficiency far surpasses the full 
BW method and is comparable to the reduced BW method. Except for periods 6 
and 10, the ST method locates the same or larger number of orbits.^ 

Unlike the double rotor map, the CHM possesses very few periodic orbits for 
small p, particularly for odd values of p. Therefore, we found that the direct appli- 
cation of the detection strategy outlined at the beginning of §3.4 would not allow 
us to complete the detection of even period orbits. Therefore, for even periods p 
we also used p + 2 as seeds and, in case of period 12, a few remaining orbits were 
located with seeds of period 15. We did not attempt to locate a maximum possible 
number of UPOs for p > 12. The numbers of such orbits (labeled with asterisks) 

^The precise reason for the failure of the ST method to detect ah period 6 and 10 orbits needs 
further investigation. We beheve that the orbits that were not detected have uncharacteristicaUy 
small convergence basins with any of the applied stabilising transformations. 
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are listed in Table 3.2 for completeness. 

As with the double rotor map, we used the symmetry of the CHM to test the 
completeness of the detected sets of orbits. It is clear from the definition of the 
CHM that all its UPOs are related by the permutation symmetry (i.e., six permu- 
tations of indices j). The column labeled "Max" in Table 3.2 lists the maximum 
number of UPOs that we were able to find using all three methods and applying the 
permutation symmetry to find any UPOs that might have been missed. As can be 
seen in Table 3.2, only a few orbits remained undetected by the ST method. 

Concluding this Section, we would like to point out that the high efficiency of 
the proposed method is primarily due to the fact that each ST constructed based 
on the stability properties of the seed orbit substantially increases the basins of 
convergence of orbits stabilised by this transformation. This is apparent in a typical 
increase of the fraction of converged seeds with the increasing value of parameter (3 
in Eq. (3.8). For example, when detecting period-10 orbits of CHM using period-12 
orbits as seeds, the fraction of seeds that converge to periodic orbits grows from 
25-30% for small P (essentially the Newton-Raphson method) to about 70% for the 
optimal value of /3. 

3.5 Summary 

In this chapter we have presented a new scheme for constructing stabilising trans- 
formations [14] which can be used to locate periodic orbits in chaotic maps with 
the iterative scheme given by Eq. (3.8). The scheme is based on the understanding 
of the relationship between the STs and the properties of eigenvalues and eigenvec- 
tors of the stability matrices of the periodic orbits. Of particular significance is the 
observation that only the unstable eigenvalues are important for determining the 
STs. Therefore, unlike the original set of transformations proposed by Schmelcher 
and Diakonos, which grows with the system size as 2"'n!, our set has cardinality of 
at most 2^, where k is the maximum number of unstable eigenvalues (i.e. the max- 
imum dimension of the unstable manifold). It is also apparent that, while the SD 
set contains a large fraction of transformations that do not stabilise any UPOs of a 
given system, all of our transformations stabilise a significant subset of UPOs. The 
dependence of the number of transformations on the dimensionality of the unstable 
manifold rather than on the system dimensionality is especially important in cases 
when we study low-dimensional chaotic dynamics embedded in a high- dimensional 
phase space. This is often the case in systems obtained from time-space discreti- 
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sation of nonlinear partial differential equations (e.g. the Kuramoto-Sivashinsky 
equation). Application of the STs approach to such high- dimensional chaotic sys- 
tems will be the subject of Chapter 4. 



Chapter IV 



Extended systems: 
kuramoto-slvashinsky equation 

The mind uses its faculty for creativity only when experience forces it 
to. 

H. J. Poincare 

In this chapter we extend the ideas presented in Chapter 3 so as to efficiently com- 
pute unstable periodic orbits (UPO) in large-scale dynamical systems arising from 
the spatial discretisation of parabolic PDEs [15]. Following the approach often 
adopted in subspace iteration methods (see [62, 88] and references therein) we con- 
struct a decomposition of the tangent space into unstable and stable orthogonal 
subspaces. On the unstable subspace we apply the method of stabilising transfor- 
mations (ST), whilst Picard iteration is performed on the complement. The method 
is extremely effective when the dimension of the unstable subspace is small com- 
pared to the system dimension. We apply the new scheme to the model example of 
a chaotic spatially extended system - the Kuramoto-Sivashinsky equation. 

4.1 Subspace decomposition 

Consider the solution of the nonlinear system 

f{x)-x = 0, xeW, f:W-^W, (4.1) 

where f{x) is assumed twice differentiable in the neighbourhood of x*, an isolated 
root of Eq. (4.1). We can approximate the solution of (4.1) by a recursive fixed point 
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procedure of the form 

x,+i = /(xi), 2 = 1,2,3,.... (4.2) 

It is well known that the iteration (4.2) converges locally in the neighbourhood of 
a solution x*, as long as all the eigenvalues of the Jacobian Df{x*) lie within the 
unit disc {z & C : \z\ < 1}. In contrast, (4.2) typically diverges if Df{x*) has an 
eigenvalue outside the unit disc. In that popular alternative - as discussed 

in Chapter 2 - is to employ Newton's method 

{Df{x,)-In)6x, = -(/(x,)-x,), (4.3) 
Xi+i = Xi + 6xi, z = l,2,3, (4.4) 

The idea of subspace iterations is to exploit the fact that the divergence of the 
fixed point iteration (4.2) is due to a small number of eigenvalues, n„, lying outside 
the unit disc. By decomposing the space M" into the direct sum of the unstable 
subspace spanned by the eigenvectors of Df{x*) 

F = Span{efc G : Df{x*)ek = XkCk, \Xk\ > 1} (4.5) 

and its orthogonal complement, Q, a modified iterative scheme is obtained. The 
application of Newton's method to the subspace P whilst continuing to use the 
relatively cheap fixed point iteration on the subspace Q, results in a highly efficient 
scheme provided dim(P) ^ dim(Q). 

To this end, let Vp G M"^"" be a basis for the subspace P C M" spanned by the 
eigenvectors of DF{x*) corresponding to those eigenvalues lying outside the unit 
disc, and Vq G R"^"'* a basis for Q, where Uu + rig = n. Then, we can define 
orthogonal projectors P and Q onto the respective subspaces, P, Q, as follows 



P = VpVp\ (4.6) 

Q = VqVj = Ir.-P (4.7) 

Note that any x G M" admits the following unique decomposition 

X = Vpp +Vqq = p + q, p:=Vpp = Px, q:=Vqq = Qx, (4.8) 



with p G M"" and q G ]R"^ Substituting (4.8) in Eq. (4.3) and multiplying the result 
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by [Vg, Vp]'^ on the left, one obtains 



V^DfVg - 




Aq 




' Vjf-q' 


V^DfV, VjDfVp - _ 




Ap 




_ vjf-p _ 



(4.9) 



Here we have used the fact that V^Vg = 0„„xns5 ^J^p = 0„^xn„, and VjDfVp = 
xn„ the latter holding due to the invariance of Df on the subspace P. Now, the 
first Hs equations in (4.9) may be solved using the following fixed point iteration 
scheme 

Ag[°l = 0, 

AgW = VjDfVgAq^^-'^ + Vjf-q, (4.10) 
Aq = AgW = 5^(V/D/V,)^(V7/-g), 

where / denotes the number of fixed point iterations taken per Newton-Raphson step. 
Since r„[V^DfVq] < 1 by construction, the iteration (4.10) will be locally convergent 
on Q in some neigbourhood of Aq - here r„[-] denotes the spectral radius. In order 
to determine Ap one solves 

{VjDfVp - 4J Ap = -Vjf + p- VjDfVgAq. (4.11) 

Note that in practice only one iteration of (4.10) is performed [62], i.e. / = 1, this 
leads to the following simplified system to solve for the correction [Ag, Ap]""" 







Aq 




' Vjf-q' 


VjDfVg VjDfVp - J„„ _ 




Ap 




_ vjf-p _ 



(4.12) 



Key to the success of the above algorithm is the accurate approximation of 
the eigenspace corresponding to the unstable modes. In order to construct the 
projectors P, Q, the Schur decomposition is used. However, primary concern of 
the work in [62, 88] is the continuation of branches of periodic orbits, where it is 
assumed that a reasonable approximation to a UPO is known. Since we have no 
knowledge a priori of the orbits whereabouts we shall need to accommodate this 
into our extension of the method to detecting UPOs. 
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4.1.1 Stabilising transformations 

An alternative approach is supplied by the method of STs, where in order to detect 
equilibrium solutions of Eq. (4.1) we introduce the associated flow 

Here g{x) = f{x)— x. With this setup we are able to stabihse all UPOs x* of 
Eq. (4.1) such that all the eigenvalues of the Jacobian Df{x*) have real part smaller 
than one. In order to stabilise all possible UPOs we study the following flow 

^ = Cg{x), (4.14) 

where C G M"^" is a constant matrix introduced in order to stabilise UPOs with 
the Jacobians that have eigenvalues with real parts greater than one. 

Substituting (4.8) in Eq. (4.13) and multiplying the result by [K/, V^]""" on the 
left, one obtains 



dq 

ds 
dp 



Thus we have replaced the original Eq. (4.13) by a pair of coupled equations, 
Eq. (4.15) of dimension Us and Eq. (4.16) of dimension n^. Since 
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and rfj\yj DgV^ < by construction, it follows that in order to detect all UPOs of 
Eq. (4.1), it is sufficient to solve 

^ = CVjg, (4.18) 

where C G M""^"" is a constant matrix. 

In [62, 88] the Schur decomposition (SchD) is used in order to construct the 
projectors P and Q. This is fine for continuation problems since one may assume 
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(a) Schur decomposition (b) Singular value decomposition 

Figure 4.1: Basins of attraction for the period— 3 orbits of the Ikeda map with 
parameter values a = 1.0, b = 0.9, k = 0.4 and rj = 6.0. Here we have chosen 
C = —1 since in this example the unstable subspace is one-dimensional. 

from the offset that they possess an initial condition Xq sufficiently close to a UPO 
such that the SchD of Df{xo) gives a good approximation to the eigenspace of 
Df{x*). However, it is well known that the eigenvectors of the perturbed Jacobian 
Df{x* + 5x) behave erratically as we increase 5x. In order to enlarge the basins 
of attraction for the UPOs we propose that singular value decomposition (SVD) be 
used instead. That is we choose an initial condition xq and construct the SVD of its 
preimage, i.e. Df{f~^{xo)) = USW'^ (or in the continuous case D(j)^{(f)~'^{xo)) = 
USW'^ for some time T), the columns of U give the stretching directions of the map 
at xq, whilst the singular values determine whether the directions are expanding or 
contracting. It is these directions which we use to construct the projectors P and 
Q. Due to the robustness of the SVD we expect the respective basins of attraction 
to increase. 

It is not necessary in practice to decompose Eq. (4.13) in order to apply the new 
ST. Rather we can express C in terms of C and Vp. To see this we add Vg times 
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Eq. (4.17) to Vp times Eq. (4.18) to get 



X] 



(4.19) 



where the second hne follows from Eq. (4.7). From this we see that the following 
choice of C is equivalent to the preceding decomposition 



C = h + Vp{C-I^^)VJ. 



(4.20) 



Thus in practice we compute Vp and C at the seed xq in order to construct C and 
then proceed to solve Eq. (4.14). 

The advantage of using the SVD rather than the SchD can be illustrated by the 
following example. Consider the Ikeda map [50]: 



/(x) := 








yi+i_ 





a + b{xi cos {(pi) - Ui sin {(pi)) 
b{xi sin(0i) + Hi cos {(pi)) 



(4.21) 



where (pi = k — ri/{l + xf + yf) and the parameters are chosen such that the map 
has a chaotic attractor: a = 1.0, b = 0.9, k = 0.4 and t] = 6.0. For this choice 
of parameters the Ikeda map possesses eight period— 3 orbit points (two period-3 
orbits and two fixed points, one of which is on the attractor basin boundary). In 
our experiments we have covered the attractor for Eq. (4.21) with a grid of initial 
seeds and solved the associated flow for p = 3, i.e., g{x) = f^{x) — x. This is done 
twice, firstly in the case where the projections P and Q are constructed via the 
SchD and secondly when they are constructed through the SVD. Since all UPOs of 
the Ikeda map are of saddle type, the unstable subspace is one-dimensional and we 
need only two transformations: C = 1 and C = —1. Figure 4.1 shows the respective 
basins of attraction for the two experiments with C = —1. It can be clearly seen 
that the use of SVD corresponds to a significant increase in basin size compared to 
the SchD. Note that with C* = — 1 we stabilise four out of eight fixed points of 
The other four are stabilised with C = 1. The corresponding basins of attraction 
are shown in Figure 4.2. Note that for this choice of ST the choice of basis vectors 
is not important, since Eq. (4.20) yields C = J, so that the associated flow is given 
by Eq. (4.13). 



4.2 Implementation 



65 




Figure 4.2: The basins of attraction for the Ikeda map for the choice of C = 1. 
Fixed points of with negative unstable eigenvalues are stable stationary solutions 
of the associated flow, while those with positive eigenvalues are saddles located at 
the basin boundaries. 



4.2 Implementation 

We now wish to apply these ideas to parabolic PDEs of the form (1-14). The 
numerical solution of such PDEs are typically based upon their representation in 
terms of a truncated system of nonlinear ODEs. Thus, we will be concerned with 
the detection of UPOs for large systems of ODEs: 

dx 

— = v{x), X G M'', (4.22) 

where v is derived from an appropriate discretisation procedure (i.e. finite differ- 
ences, finite element, spectral method) and n ^ 1. 

A typical approach in the determination of UPOs for flows is via a Poincare 
surface of section (PSS). By "clever" placement of an {n — l)-dimensional manifold 
in the phase space, the problem is reduced to a discrete map defined via intersections 
with the manifold. However, a correct choice of PSS is a challenging problem in 
itself. Due to the complex topology of a high-dimensional phase space, the successful 
detection of UPOs will be highly dependent upon the choice of section. When the 
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choice of a suitable PSS is not obvious a priori, we found it preferable to work with 
the full flow, adding an auxiliary equation to determine the integration time T. 

Let X denote the flow map of Eq. (4.22). Then we define the associated 

flow as follows 



— = Cg 
ds 



(4.23) 



The additional equation for T is constructed such that T is always changing in the 
direction that decreases Id'^ix) — x\, i.e. 



dT 
ds 



oc 



d\gl 
dT ' 



or, more precisely 



dT 
ds 



—av[ 



\x))-{<f{x)-x). 



(4.24) 



d 


X 




C(0^(x) - x) 


ds 


T 




—av{(j)'^{x)) ■ {4>'^{x) — x) 



Here a > is a constant which controls the relative speed of convergence of 
Eq. (4.24). This leads to the following augmented flow which must be solved to 
detect UPOs of (4.22): 



(4.25) 



Note that the augmented flow (4.25) is derived by integrating a nonlinear PDE - in 
our case the KSE - for some time T and will become increasingly stiff for larger T. 

Several approaches have been proposed for the solution of stiff systems of ODEs; 
see for example, the review by Shampine and Gear [87]. Of all these techniques the 
general-purpose codes contained within the ODEPACK software package [44] are re- 
garded as some of the best available routines for the solution of such systems. Thus, 
in our numerical experiments we use the stiff solver dlsodar from the ODEPACK 
toolbox to integrate (4.25). dlsodar is a variable step-size solver which automat- 
ically changes between stiff and nonstiff modes. In particular, as we approach a 
steady state of Eq. (4.25) dlsodar will take increasingly larger time-steps, leading 
to super linear convergence in the neighbourhood of the solution. 

To use the solver dlsodar, we must provide a routine that returns the value 
of the vector field (4.25) evaluated at a given point (x,T). Here we need the so- 
lution of Eq. (4.22) which is obtained by applying a suitable numerical integration 
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scheme; see the next section and Appendix E for further details. The ODEPACK 
software package makes use of the Jacobian matrix of the system being integrated 
and provides the option of computing the Jacobian via finite differences or via a user 
supphed routine. Note that in the case that the flow is expected to be stiff much of 
the time, it is recommended that a routine for the Jacobian is supplied and we do 
this. The derivative of (4.25) with respect to {x,T) is given by 

C{Jt - /„) CvT 
—a{v]^{JT — In) + 9^ DvJt) —a{g^ Dvvt + v^vt) 

where Jt = d(j)^{x)/dx, vt = v{(j)'^{x)), Dv = dv/dx and g = (j)^{x) — x as usual. 

Quite often one might wish to terminate simulation before the usual stopping 
criteria of, for example, a maximum number of steps taken or certain tolerances 
having been reached. A useful feature of the dlsodar algorithm is that it allows the 
optional user supplied routine to do just this. To be more exact, it determines the 
roots of any of a set of user supplied functions 

and returns the solution of (4.25) at the root if it occurs prior to the normal stopping 
criteria. 

Note that, to increase the efficiency of the algorithm we wish to avoid the fol- 
lowing two instances: firstly, due to the local nature of the STs, we should stop the 
search if we wander too far from the initial condition, and secondly, since our search 
is governed by the dynamics of Eq. (4.23) and not by those of (4.22), we might 
move off the attractor after some time period so the convergence to a UPO becomes 
highly unlikely. In our numerical experiments we supply the following function 

hi = a-\g\, (4.27) 

where a G M is a constant and | ■ | denotes the L2 norm. In practice, we have found 
that there exists a threshold value of a, such that convergence is highly unlikely 
once the norm of g surpasses it. Note that we also restrict the maximum number of 
allowed integration steps since the convergence becomes less likely if the associated 
fiow is integrated for a long time. 

We must also provide two tolerances, rtol and atol, which control the local error 
of the ODE solver. In particular, the estimated local error in X = {x, T) will be 



(4.26) 
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controlled so as to be less than 

rtol ■ I |X| loo + atol. 

Thus the local error test passes if, in each component, either the absolute error is 
less than atol or the relative error is less than rtol. The accuracy with which we 
would like to solve the flow (4.25) is determined by the stability properties of (4.22). 
To understand this, we note that in evaluating the RHS of (4.25) it is the solution of 
Eq. (4.22) at time T, i.e. (f)^{x), which is critical for error considerations. Suppose 
that our initial point lies within 6x of a true trajectory x. Then (f)'^{x + 6x) lies 
approximately within e^'^6x of the true trajectory, (fP^^x), where A is the largest 
Lyapunov exponent of the system. Since A is positive for chaotic systems, the error 
grows exponentially with the period and we should take this into account when 
setting the tolerances rtol and atol. This leads us to the following settings for the 
tolerances 

rtol = atol = 10-Ve^^°, (4.28) 

where Tq is the initial period and A is the largest Lyapunov exponent of the flow v. 
We have computed the Lyapunov exponent using the algorithm due to Benettin et 
al [2]; see Appendix F for a description of the routine as well as a brief review of 
Lyapunov exponents. 

4.2.1 Kuramoto-Sivashinsky equation 

We have chosen the Kuramoto-Sivashinsky equation (KSE) for our numerical ex- 
periments. It is the simplest example of spatiotemporal chaos and has been studied 
in a similar context in [9, 57, 100], where the detection of many UPOs has been 
reported. We work with the KSE in the form 

ut = --{u^)x-u,.-u,xxx, (4.29) 

where x G [0, L] is the spatial coordinate, t G is the time and the subscripts x, 
t denote differentiation with respect to space and time. For L < 2n, u{x,t) = is 
the global attractor for the system and the resulting long time dynamics are trivial. 
However, for increasing L the system undergoes a sequence of bifurcations leading 
to complicated dynamics; see for example [54]. 

Our setup will be close to that found in [57]. In what follows we assume periodic 
boundary conditions: u{x,t) = u{x + L,t), and restrict our search to the subspace 
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of antisymmetric solutions, i.e. u{x,t) = —u{L — x,t). Due to the periodicity of 
the solution, we can solve Eq. (4.29) using the pseudo-spectral method [36, 94]. 
Representing the function u{x, t) in terms of its Fourier modes: 

u{x, t) := [n] = J2 ^ke-''"^'', (4.30) 

where 

-L 



U := [. . . ,U_i,Mo,Ml, 



y , Uk := T[u]k = y / u{x, t)e"'''^dx, (4.31) 

Jo 



we arrive at the following system of ODEs 

^ = [{kqr - {kqr]uk + '^H^-'[u])% ■ (4.32) 

Here q = 27r/L is the basic wave number. Since u is real, the Fourier modes are 
related by u^k = Wfc- Furthermore, since we restrict our search to the subspace of 
odd solutions, the Fourier modes are pure imaginary, i.e. y{t{uk) = 0. 

The above system is truncated as follows: the Fourier transform JF is replaced 
by its discrete equivalent 

N-l ^ N~l 

au ■■= r^VAk = J2 ^i^jy^'' ' ^i^j) ■■= ^n'K = ]^ E ".e"''""^' , (4.33) 

j=0 k=0 

where xj = L/N and ajy^k = clI- Since ao = due to Galilean invariance and 
setting a7v/2 = (assuming N is even), the number of independent variables in the 
truncated system is n = N/2 — 1. The truncated system looks as follows: 

dk = [{kqf - (A;g)>, + '-!^J^j,[(J^^^[a])% , (4.34) 

with k = 1, . . . ,n, although in the Fourier transform we need to use ak over the full 
range of k values from to — 1. 

The discrete Fourier transform jF/v can be computed using fast Fourier transform 
(FFT). In Fortran and C, the routine REALFT from Numerical Recipes [77] can be 
used. In Matlab, it is more convenient to use complex variables for a^. Note that 
Matlab function f ft is, in fact, the inverse Fourier transform. 

To derive the equation for the matrix of variations, we use the fact that is a 
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linear operator to obtain 
ddk 



[{kqf - {kq)^]6kj + tkqJ^N[J^N'H ® ^N^i^kj]] , 3 = l,...,n, (4.35) 



where ® indicates componentwise product, and the inverse Fourier transform is 
apphed separately to each column of S^j- Here, 6k j is not a standard Kronecker 
delta, but the N x n matrix: 



/ 
1 
1 



Okj 



\ 







V 



-1 ... 
-1 ...y 



(4.36) 



with index k running from to — 1. 

In practice the number of degrees of freedom n should be sufficiently large so 
that no modes important to the dynamics are truncated, whilst on the other hand, 
an increase in n corresponds to an increase in computation. To determine the order 
of the truncation in our numerical experiments, we initially chose n to be large and 
integrated a random initial seed onto the attractor. By monitoring the magnitude 
of the harmonics an integer k was determined such that aj < 10~^ for j > k. The 
value of n was then chosen to be the smallest integer such that: (i) n > k, and 
(ii) N = 2n + 2 was an integer power of two. The second condition ensures that 
the FFT is applied to vectors of size which is a power of two resulting in optimal 
performance. 

Note that in the numerical results to follow we work entirely in Fourier space and 
use the ETDRK4 time-stepping to solve (4.34) and (4.35); for further details con- 
cerning the method of exponential-time-differencing see Appendix E. In particular 
the method uses a fixed step-size {h = 0.25 in our calculations) thus it is necessary 
to use an interpolation scheme in order to integrate up to arbitrary times. In our 
work we have implemented cubic interpolation [77]. More precisely, to integrate up 
to time t G [ti,ti + h], where the ti are integer multiples of the step-size h. We 
construct the unique third order polynomial passing through the two points a(ti) 
and a{ti + h), with derivatives a'(tj) and a'{ti + h) at the respective points. In this 
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way we obtain the following cubic model: 



p{s) = [2a{ti) + a'{ti) + a{ti + h) - 2a{ti + h)]s^ + [3a{ti + h) - 3a{ti) 
- 2a'{ti) - a\ti + h)]s^ + a'{ti)s + a(ti), 



(4.37) 



where the parameter s = (t — ti)//i G [0, 1]. 
4.2.2 Numerical results 

We now present the results of our numerical experiments. We begin with a com- 
parison between our method and the nonlinear least squares solver Imder from the 
MINPACK software package [66]. Note that Imder is an implementation of the 
Levenburg-Marquardt algorithm - see §2.3. We have chosen it because it has re- 
cently been applied successfully to detect many UPOs of the closely related complex 
Ginzburg-Landau equation. 

In order to determine UPOs via the Imder routine we introduce the following 
augmented system 



where the second equation is motivated in analogous fashion to the auxiliary equa- 
tion (4.24). To use the Imder solver we must provide a routine that returns the 
value of F evaluated at a given point (a;, T). Also, we provide a routine to compute 
the Jacobian analytically rather than use finite differences, the formula of which 
differs from Eq. (4.26) only by multiplication by a constant matrix. In addition to 
this the user must supply three tolerances: ftol, xtol, and gtol. Here, ftol measures 
the relative error desired in the sum of squares, xtol the relative error desired in 
the approximate solution, and gtol measures the orthogonality desired between the 
vector function F and the columns of the Jacobian. In the computations performed 
in the next section we set the tolerances to the recommended values: 



Finally, we specify a maximum number of function evaluations allowed during each 
run of the Imder algorithm in order to increase its efficiency. 





ftol = 10 



8 



Xtol = 10"^ and gtol = 0. 
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For our method we use the set of matrices proposed in §4.1.1 with 



C = Ci 



■SD, 



(4.39) 



since within the low-dimensional unstable subspace it is possible to apply the full 
set of Schmelcher-Diakonos (SD) matrices. The UPOs determined from our search 
will then be used as seeds to determine new cycles. Here we proceed in analogous 
fashion to Chapter 3 by constructing STs from the monodromy matrix, D(j)^*{a*), 
of the cycle (a*,T*). We then solve the augmented flow (4.25) from the new initial 
condition (a*,T), where the time T is chosen such that 



Note that any given cycle may exhibit many close returns, particularly longer cycles, 
thus in general a periodic orbit may produce many new initial seeds. This is an 
especially useful feature, since we do not have to recompute the corresponding STs. 

Initially, to determine that a newly detected cycle, (x*,T*), was different from 
those already found, we first checked whether the periods differed, that is, we de- 
termined whether |T* — T'| > Ttoi for all previously detected orbits. If two orbits 
where found to have the same period, we then calculated the distance between the 
first components of x*, and all other detected orbits y*. However, in practice we 
found, that if two orbits have the same period, then either they are the same or they 
are related via symmetry, recall that if u{x,t) is a solution then so is —u{L — x,t). 
Thus, in order to avoid the convergence of UPOs that are trivially related by sym- 
metry we will consider two orbits as being equal if their periods differ by less than 
the tolerance Ttoi. Of course, this criterion can, in theory, lead us to discard cycles 
incorrectly, however, this is highly unlikely in practice. 

Stabilising transformations vs Levenberg-Marquardt 

The search for UPOs is conducted within a rectangular region containing the chaotic 
invariant set. Initial seeds are obtained by integrating a random point within the 
region for some transient time r. Once on the attractor, the search for close returns 
within chaotic dynamics is implemented. That is, we integrate the system from the 
initial point on the attractor until a(to) ~ ^(ti) for some times to < ^i, and use the 
close return, (a(to),To), where Tq = ti — t^, as our initial guess to a time-periodic 
solution. In order to build the STs we solve the variational equations for each seed 
starting from the random initial point, a(0), for time r + to- In order to construct 



a*(0) ^ a*(r), and T (modT*) ^ 0. 



(4.40) 



4.2 Implementation 



73 



Table 4.1: The number of distinct periodic solutions for the Kuramoto-Sivashinsky 
equation detected by the method of STs. Here L = 38.5 and a = 0.25. 



Period 


C 


iVpo 




iVfcv 




Work 




Co 


28 


498 


252 


10 


412 


10 - 100 




16 


296 


684 


32 


1196 






44 


794 


936 


42 


1608 




Co 


235 


395 


903 


78 


2151 


100 - 250 




64 


256 


1294 


112 


3086 






299 


651 


2197 


190 


5237 



Table 4.2: The number of distinct periodic solutions for the Kuramoto-Sivashinsky 
equation detected by the Levenberg-Marquardt algorithm Imder with L = 38.5. 



Period 


iVpo 


iVhit 


A^fev 




Work 


10 - 100 


42 


497 


20 


15 


260 


100 - 250 


208 


291 


500 


452 


7732 



the matrix Vp, we apply the SVD to the matrix 

D(j)^+''>{a{0)) = USW'^, (4.41) 

the corresponding ST is given by 

C = I^ + VpiC-I^jV;, (4.42) 

where Vp = Ujk, j = 1, . . . ,n, k = 1, . . . , riu, i.e. the first columns of U in (4.41). 
Here Uu is the number of expanding directions which is determined by the number 
of singular values with modulus greater than one. 

We examine two different system sizes: L = 38.5, for which the detected UPOs 
typically have one positive Lyapunov exponent, and L = 51.4, for which the UPOs 
have either one or two positive Lyapunov exponents. The corresponding systems 
sizes are n = 15 and n = 31 respectively. Our experiments where conducted over two 
separate ranges. We began by looking for shorter cycles with period T G [10, 100], 
the lower bound here was determined a posteriori so as to be smaller than the 
shortest detected cycle. We then searched for longer cycles, T G [100, 250] to be 
more precise, where the maximum of T = 250 was chosen in order to reduce the 
computational effort. 
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In our calculations we set the positive constant a = 0.25 in Eq. (4.25). Using the 
solver dlsodar we integrated 500 random seeds over both ranges for time s = 150, 
the seeds where chosen such that |0"^"(a(O))— a(0)| < 1.0. If the flow did not converge 
within 1000 integration steps, we found it more efficient to terminate the solver and 
to re-start with a different ST or a new seed. As mentioned in the previous section 
we choose a constant a = 50 experimentally so that integration is terminated if the 
norm of g grows to large, i.e. \g{x)\ = \(f)'^{x) — x\ > a. 

Typically on convergence of the associated flow the UPO is determined with 
accuracy of about 10^'' at which point we implement two or three iterations of the 
Newton-Armijo rule to the following system 

Mt - In V{(f)^{x)) \ ( Sx\ _ f (f)^{x) - X 

v{x) ) \5T ) ~ \ 

in order to allow convergence to a UPO to within roundoff error. The linear system 
(4.43) was introduced by Zoldi and Greenside [100] as a method for determining 
UPOs for extended systems in its own right, and in particular, has been applied 
successfully in the detection of UPOs for the KSE. 

Similarly, we run the Imder routine from the same 500 seeds over the two different 
time ranges. The routine terminates if one of the following three scenarios arise: 
(i) a predefined maximum number of function evaluations is exceeded, we set the 
maximum number of function evaluations equal to 1000, (ii) the error between two 
consecutive steps is less than xtol, but the sum of squares is greater than ftol, 
indicating a local minimum has been detected, or (iii) both xtol and ftol are satisfied 
indicating that convergence to a UPO has been obtained. 

In order to make a comparison between the efficiency of the two methods we 
introduce the measure of the work done per seed: 

Work = A/'fev + n X iVjev (4.44) 

Here A^fev is the average number of function evaluations per seed, A'jev the average 
number of Jacobian evaluations per seed and n is the size of the system being solved. 
The expression in (4.44) takes into account the fact that evaluation of the Jacobian 
is n times more expensive than evaluation of the function itself. 

We present the results of our experiments in Tables 4.1 - 4.4. Here A^po denotes 
the number of distinct orbits found, Ahit gives the number of times we converged to 
a UPO, and A^fev, Ajev, and Work, are as defined above. In Tables 4.1 and 4.3 the 
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Figure 4.3: Illustration of two UPOs of KSE detected from a single seed. We 
show both a level plot for the solutions and a projection onto the first two Fourier 
components. Since u{x,t) is antisymmetric on [0,L], it is sufficient to display the 
space-time evolution of u{x, t) on the interval [0, L/2]: (a) Seed with time T = 37.0, 
(b) a periodic solution of length T = 36.9266 detected with stabilising transfor- 
mation C = +1 and (c) a periodic solution of length T = 25.8489 detected with 
stabilising transformation C = —1. 



performance of the stabilising transformations are analysed both collectively and on 
an individual basis; here the different Cj denote the different SD matrices. In Table 
4.1 Co = +1 and Ci = —1, whilst in Table 4.3 we have 
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-1 


-1 



The total work done per seed is given by the sum over all the SD matrices and is 
denoted by {Cj}. 

For L = 38.5 we detected a total of 343 UPOs using the ST method as compared 
to 250 UPOs using the Imder algorithm, whilst for L = 51.4 we found 144 and 
106 UPOs with the respective methods. Both methods found roughly the same 
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Table 4.3: The number of distinct periodic solutions for the Kuramoto-Sivashinsky 
equation detected by the method of stabilising transformations. Here L = 51.4 and 
a = 0.25. 
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n 
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1 71 




00 
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5 


174 


666 


36 


1818 






2 


139 


496 


36 


1648 






30 


1417 


3885 


205 


10445 




Co 


51 


330 


628 


17 


1172 






7 


209 


807 


27 


1671 




C2 


17 


138 


936 


31 


1928 






21 


177 


917 


26 


1749 


100 - 250 


C4 


11 


161 


877 


36 


2029 




C5 





116 


960 


43 


2336 




c. 


1 


117 


975 


42 


2319 






6 


161 


879 


35 


1999 






114 


1409 


6979 


257 


15203 



number of orbits when searching for shorter period cycles. However, as can be seen 
by comparing the work done per seed in Tables 4.1 - 4.4 Imder performs better 
than the ST method in that case. The situation changes, however, when we look 
at the detection of longer cycles. For L = 38.5 in particular, we see that the ST 
method computes many more orbits with a considerable reduction in the cost. In 
fact, using the identity matrix alone, the ST method detects more UPOs than Imder 
at approximately a quarter of the cost. For larger system size, L = 51.4, the ST 
method still detects more orbits than Imder. However, in doing so a considerable 
amount of extra work is done. It is important to note here, that the increase in work 
is not due to any deficiency in the ST method. Rather, the rise in the computational 
cost is brought about due to the fact that not all the SD matrices work well. For 
example, the subset of matrices {Co, C2, C3, C4} C Csd, detect approximately 90% 
of the longer period UPOs, yet they account for less than 50% of the overall work. 
Indeed, for this particular choice of matrices the ST method still detects more orbits, 
but more importantly, its performance now exceeds that of Imder. 
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Table 4.4: The number of distinct periodic solutions for the Kuramoto-Sivashinsky 
equation detected by the Levenberg-Marquardt algorithm Imder with L = 51.4. 



Period 


iVpo 


iVhit 




N- 


Work 


10 - 100 


20 


480 


64 


52 


1728 


100 - 250 


86 


356 


374 


332 


10998 



Note that it is not surprising that the SD matrices do not perform equally well, 
this is, after all, what we would have expected based upon our experience with maps. 
However, it is still important to try all SD matrices - when possible - in order to 
compare their efficiency, especially if one wishes to construct minimal sets of ST 
matrices. 

For detecting longer period orbits the performance of Imder suffers due to the 
large increase in the number of Jacobian evaluations needed; see Tables 4.2 and 4.4. 
Recall, that Imder is an implementation of the Levenberg-Marquardt algorithm, 
and that in particular, it computes the Jacobian once on each iterate. Thus, we 
may view the number of Jacobian evaluations as the number of steps required in 
order to converge. The increase in the number of iterations necessary to obtain 
convergence can be understood by considering the dependance of \(p'^{x) — x\ on 
X. For increasing period the level curves of l^fp are squeezed along the unstable 
manifold of (f)'^{x), resulting in a complicated surface with many minima, both local 
and global, embedded within long, winding, narrow "troughs". Note that this is 
a common problem with all methods that use a cost function to obtain "global" 
convergence, since these methods only move in the direction that decreases the cost 
function. 

This can be explained by the following heuristics: for simplicity let us assume 
that we are dealing with a map Xn+i = f{xn), whose unstable manifold is a one- 
dimensional object. In that case, we may define a one- dimensional map, locally, 
about a period-p orbit, x*, of the map / as 

h{s) = \g{x* + Sx)\'^. (4.45) 

Here g = f^{x) — x as usual, 5x = x{s) — x* is small, and we only allow x{s) to 
vary along the unstable manifold. Expanding g in a Taylor series about the periodic 
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orbit, x*, we obtain 

h{s) = \g{x*) + Dg{x*)5x + 0{Sx^)\'^ 

= \iDr{x*)-In)6x + Oi5x'')\^ 

= 6x'^{AP - Inf6x + 0{6x^), 

where the third hne follows since 6x is an eigenvector of Df{x*) with corresponding 
eigenvalue A, and A = diag(A, ... A). Note that in the above we assume that p is 
large but finite so that the term (A^ — remains bounded. Now, close to the 
periodic orbit, h is approximately a quadratic form with slope of the order A^, and 
since |A| > 1, it follows that if we move along the unstable manifold {gl"^ will grow 
quicker for larger periods, or, in other words, the level curves are compressed along 
the unstable direction. 

Since the Imder algorithm reduces the norm of / at each step, it will typically 
follow the gradient to the bottom of the nearest trough, where it will start to move 
along the narrow base towards a minimum. Once at the base of a trough, however, 
Imder is forced into taking very small steps, this follows due to the nature of the 
troughs, i.e. the base is extremely narrow and winding, and since Imder chooses 
its next step-size based on the straight line search. One of the benefits of the ST 
method is that it does not need to decrease the norm at each step and thus, does 
not suffer from such considerations. 

Another advantage of our method is that we can converge to several different 
UPOs from just one seed, depending upon which ST is used. Figure 4.3 shows one 
of such cases, where Eq. (4.25) is solved from the same seed for each of the 2"" 
STs {uu = 1 in this example), with two of them converging to two different UPOs. 
Figure 4.3a shows the level plot of the initial condition and a projection onto the first 
two Fourier components. Figures 4.3b and 4.3c show two unstable spatiotemporally 
periodic solutions which where detected from this initial condition, the first of period 
T = 36.9266 was detected using C = +1, whilst the second of period T = 25.8489 
was detected using C = — 1. The ability to detect several orbits from one seed 
increases the efficiency of the algorithm. 

Seeding with UPOs 

In order to construct a seed from an already detected orbit, (a*,T*), we begin by 
searching for close returns. That is, starting from a point on the orbit, a* (to), we 
search for a time ti such that a{ti) is close to a(to)- As long as TmodT 7^ 0, where 
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Table 4.5: Number of distinct orbits detected using the method of stabihsing trans- 
formations with periodic orbits as seeds The number of seeds is 489. L = 38.5, 
a = 0.25. 



c 


iVpo 


iVhit 


A^fev 




Work 


Co 


46 


209 


2814 


136 


4854 




52 


198 


2819 


130 


4769 




98 


407 


5633 


266 


9623 



Table 4.6: Number of distinct orbits detected using the method of stabilising trans- 
formations with periodic orbits as seeds. The number of seeds is 123. L = 51.4, 
a = 0.25. 
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iVpo 


iVhit 






Work 


Co 


1 


21 


2797 


48 


4333 


Ci 


2 


37 


2815 


40 


4095 


C2 








216 


3 


312 







1 


187 


2 


251 




3 


59 


6015 


93 


8691 



T = ti — to, we take {a*(to),T) as our initial guess to a time periodic solution. For 
longer period cycles, we can find many close returns by integrating just once over 
the period of the orbit. Shorter cycles, however, produce fewer recurrences and, 
in general, must be integrated over longer times to find good initial seeds. In our 
experiments we searched for close returns T G (10, 2T*) if T* < 150.0; otherwise, 
we chose T G (10, T*). 

Stabilising transformations are constructed by applying the polar decomposition 
to the matrix G = QB, where G is defined by 

G := V{SA-ln)V-\ 

Here V and A are defined through the eigen decomposition of Dcp'^* (a*{to)) = 
VAV~^, and S = diag(±l, ±1, . . . , ±1). The different transformations are then 
given by C = — Q""". Recall from Chapter 3 that changing the signs of the stable 
eigenvalues is not expected to result in a substantially different stabilising transfor- 
mation, thus we use the subset of S such that Sa = 1 for i > n„. For L = 38.5, we 
have just two such transformations, since all periodic orbits have only one unstable 
eigenvalue. Whilst for L = 51.4, we will have either two or four transformations 
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depending upon the number of unstable eigenvalues of Dcp"^* {a*{to)). 

In the calculations to follow, a close return was accepted if \a*(ti) — a*(to)| < 2.5. 
Note that, this was the smallest value to produce a sufficient number of recurrences 
to initiate the search. A particular cycle may exhibit many close returns. In the 
following, we set the maximum number of seeds per orbit equal to 5. To obtain 
convergence within tolerance 10~^ we had to increase the integration time to s = 250 
and the maximum number of integration steps allowed to 2000. If we converged to 
a UPO then as in the preceding section, we apply two or three Newton- Armijo steps 
to the linear system (4.43) in order to converge to machine precision. 

The results of our experiments are given in Tables 4.5 and 4.6. As in the previous 
section, A^po denotes the number of distinct orbits found, Ahit the number of times 
we converged to a UPO, A^fev the average number of function evaluations per seed, 
and Ajac the average number of jacobian evaluations per seed. The computational 
cost per seed is measured in terms of the average number of function evaluations 
per seed and is defined as in Eq. (4.44). In Tables 4.5 and 4.6 the different Ci can 
be uniquely identified by the signature of pluses and minuses defined through the 
corresponding matrix 5*. For L = 38.5, the matrices Co and Ci correspond to the 
signatures (+, ...,+) and (— , +,...,+), respectively, whilst for L = 51.4, the ma- 
trices Co, Ci, C2 and C3 correspond to (+,..., +), (-,+,..., +), (-, -,+,..., +), 
and (+, —,...,-!-) respectively. As before, the total work done per seed is defined 
as the sum over all matrices and is denoted by {Ci}. 

For L = 38.5 we where able to construct 489 seeds from 343 previously detected 
periodic orbits. From the new seeds we found a further 98 distinct UPOs, bringing 
the total number of distinct orbits detected for L = 38.5 to 433. An important 
observation to make is that both matrices perform equally well, as can be seen from 
Table 4.5. This is in contrast to the SD matrices for which the performance of the 
different matrices varies greatly; see Tables 4.1 and 4.3. The result of using periodic 
orbits as seeds has, however, increased the cost per seed as compared with the cost 
of seeding with close returns within a chaotic orbit, see Tables 4.1 and 4.5. The 
increased computation is due mainly to the fact that the close returns obtained 
from periodic orbits are not as good as those obtained from a chaotic orbit. Recall 
that the stabilising transformations are based on the local invariant directions of the 
orbit, and we would expect the performance to suffer as we move further from the 
seed. 

For larger system size, where the system becomes more chaotic, the construc- 
tion of seeds from close returns becomes increasingly difficult. For L = 51.4, we 
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Figure 4.4: Illustration of how a UPO can be used as a seed to detect new cycles. 
We show the projection onto the first two Fourier components of the initial seed 
(T = 237.4470) and the three detected orbits. Here L = 38.5 and n = 15. 



constructed 123 initial seeds from the 144 periodic orbits using the method of near 
recurrences. Here we detected only 3 new distinct UPOs. It is important to note, 
however, that although we do not find many new UPOs for the case L = 51.4, the 
method is still converging for approximately 20% of all initial seeds. Also, in Table 
4.6, the poor performance of the matrices C2, C3 as compared to that of Co, Ci, is 
due to the fact that only a small percentage of seeds were constructed from orbits 
with two positive Lyapunov exponents. 

One advantage of using periodic orbits as seeds is that we can construct many 
seeds from a single orbit. Figure 4.4 gives an illustration of three periodic orbits 
which were detected from a long periodic orbit used as seed. Figure 4.4 shows 
the projection onto the first two Fourier components of the initial seed (of period 
T = 237.4470) and the three orbits detected whose periods in descending order 
are T = 200.428, T = 58.7515 and T = 52.1671. For each seed determined from 
a particular periodic orbit the stability transformations are the same. Hence, the 
ability to construct several seeds from a single orbit, increases the efficiency of the 
scheme. 
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4.3 Summary 

In this chapter we have presented a scheme for detecting unstable periodic orbits 
(UPOs) in high-dimensional chaotic systems [13, 15] based upon the stabilising 
transformations (STs) proposed in [19, 83]. Due to the fact that high- dimensional 
systems studied in dynamical systems typically consist of low-dimensional dynam- 
ics embedded within a high-dimensional phase space, it is possible to increase the 
efficiency of the STs approach by restricting the construction of such transforma- 
tions only to the low- dimensional unstable subspace. Following the approach often 
adopted in subspace iteration methods [62], we construct a decomposition of the 
tangent space into unstable and stable orthogonal subspaces. We find that the use 
of singular value decomposition to approximate the appropriate subspaces is prefer- 
able to that of Schur decomposition, which is usually employed within the subspace 
iteration approach. As illustrated with the example of the Ikeda map, the decompo- 
sition based on singular value decomposition is less susceptible to variations in the 
properties of the tangent space away from a seed and thus produce larger basins of 
attraction for stabilised periodic orbits. Within the low-dimensional unstable sub- 
space, the number of useful STs is relatively small, so it is possible to apply the full 
set of Schmelcher-Diakonos matrices. The detected orbits were then used as seeds in 
order to search for new UPOs, thus enabling us to apply the STs introduced in [14]. 



Chapter V 
Summary and outlook 



In mathematics you don't understand things. You just get used to them. 
J. von Neumann 

5.1 Summary 

In this thesis we have discussed in detail the appUcation to high- dimensional systems 
of the method of stabilising transformations (ST). The main advantage of the ST 
approach as compared to other methods is its excellent convergence properties. In 
particular, the basins of attraction are much larger than the basins produced by other 
iterative schemes, and are simply connected regions in phase space. The convergence 
properties for methods of Newton or secant type, can be improved upon by using 
a line search, the step-size being determined by a suitable cost function. However, 
such methods have no way of differentiating between true roots and local minima 
of the cost function; a problem which increases significantly with system dimension 
due to the topology of multi-dimensional flows. 

The application to high-dimensional systems however, is not straightforward. 
The set of transformations conjectured by Schmelcher and Diakonos (SD) have two 
major failings: (i) the complete set of transformations has cardinality 2"n! (here n 
is the system dimension), and (ii) this set contains a certain redundancy, in that, 
not all matrices are useful. The main focus has been on trying to understand the 
theoretical foundations of the stabilising transformations in order to determine a 
new minimal set of such transformations enabling efficient detection for systems 
with n> A. 

Our approach is based on the understanding of the relationship between the sta- 
bilising transformations and the properties of the eigenvalues and eigenvectors of the 
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stability matrices of the periodic orbits. Of particular significance is the observation 
that only the unstable eigenvalues are important for determining the stabilisation 
matrices. Therefore, unlike the SD matrices, whose numbers grow with system di- 
mensionality as 2"n!, our set has at most 2^^ elements, where k is the dimension of 
the unstable manifold. The dependence of the number of transformations on the 
dimensionality of the unstable manifold is especially important in cases when we 
study low-dimensional chaotic dynamics embedded within a high- dimensional phase 
space. 

The ST method was originally developed for the detection of periodic orbits 
in time-discrete maps. However, the periodic orbits can be used similarly to infer 
the properties of time-continuous dynamical systems. Here, the ST approach is 
typically applied directly to the corresponding Poincare map, where the Poincare 
surface of section (PSS) can be a stroboscopic map or defined in terms of phase space 
intersections. For high- dimensional systems, however, the correct choice of a PSS is 
a challenging problem in itself. Due to the complex topology of a high- dimensional 
phase space, the successful detection of unstable periodic orbits (UPOs) will depend 
upon the choice of PSS. A major advantage of our method is that it does not require 
a PSS in order to apply the ST method to determine UPOs for flows. 

One apparent drawback of the new scheme is that a small set of UPOs needs 
to be available for the construction of the stabilising transformation at the start 
of the detection process. However, we have shown that it is possible to construct 
the stabilising transformations without the knowledge of UPOs. Recalling that 
the stabilising transformations depend mostly on the properties of the unstable 
subspace, we use the fact that the decomposition into stable and unstable subspaces 
can be defined at any, not just periodic, point on the chaotic set. The decomposition 
is done in a process similar to that used in the subspace iteration algorithm [62], 
the full set of Csd matrices can then be applied in the low- dimensional unstable 
subspace. 

The new transformations were tested on a kicked double rotor map, three sym- 
metrically coupled Henon maps and the Kuramoto-Sivashinsky equation (KSE). For 
the time-discrete maps our aim was to achieve the detection of plausibly complete 
sets of periodic orbits of low periods up to as high a period as was computationally 
feasible. In both cases our algorithm was able to detect large numbers of UPOs 
with high degree of certainty that the sets of UPOs for each period were complete. 
We have used the symmetry of the systems in order to test the completeness of the 
detected sets. On the other hand, when the aim is to detect as many UPOs as 



5.2 Outlook 



85 



possible without verifying the completeness, the symmetry of the system could be 
used to increase the efficiency of the detection of UPOs: once an orbit is detected, 
additional orbits can be located by applying the symmetry transformations. 

In the case of the KSE our goal was somewhat different [13, 15]. Here the aim 
was to construct and implement stabilising transformations determined from an 
arbitrary point in phase space. The proposed method for detecting UPOs has been 
tested on both a 15 and 31- dimensional system of ODEs representing odd solutions 
of the KSE. Using the new set of stabilising transformations we have been able 
to detect many periodic solutions of the KSE in an efficient manner. The newly 
detected orbits were then used as seeds to detect new orbits by extending the ideas 
in [14] to the continuous case. 

5.2 Outlook 

The method of stability transformations is a powerful tool in the detection of unsta- 
ble periodic orbits for a large class of dynamical systems. Until now, its application 
has been restricted to low-dimensional dynamical systems. In this thesis we have 
used the insights gained from our two-dimensional analysis in order to extend the 
ST method efficiently in higher dimensions. However, the understanding of the the- 
oretical foundations in more than two dimensions remains a challenge. For example, 
the SD Conjecture 1.1 still remains an open problem for n > 2. One possible ap- 
proach to this question would be to show that for any orthogonal matrix Q there 
exists at least one matrix, C G Csd, which is close to Q for all n. Here the measure 
of closeness is that defined in [14]. 

The detection of spatial and temporal patterns in the time evolution of nonlinear 
PDEs is an area of research currently receiving a lot of attention; see [52, 97, 95] 
and references therein. Here, an important question is whether or not the results of 
the periodic orbit theory generalise, and more importantly, remain useful for such 
high- dimensional dynamical systems. The numerical treatment of such systems, 
however, is limited due to computational constraints. Thus, the construction of 
efficient tools for detecting spatiotemporal patterns is extremely important. Our 
numerical results for the KSE are a first step in this direction. In the future, an 
exploration of the full solution space should be conducted; the restriction to odd 
subspace is a computational convenience, rather than a physically meaningful con- 
straint. Importantly, for the full system, solutions are invariant under translations 
along the x-axis, i.e. if u{x,t) is a solution then so is u{x + A,t), VA G M. Now 
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UPOs only tell part of the story and it is expected that relative periodic solutions, 
i.e. u{x + A,t + T) = u{x,T), will become important for a proper description of 
the dynamics. One way to eliminate this arbitrariness is to supply an extra equa- 
tion, this may be done, for example, in analogy to the auxiliary equation (4.24) of 
Chapter 4: 

where g{x, A, T) = + A) — x. 

Also, it is clear that as we start to study more and more complex systems, that 
it becomes increasingly difficult to find UPOs from which to initialise the detec- 
tion process. Thus, the idea of constructing matrices from an arbitrary point in 
phase space becomes increasingly important, and although numerical results indi- 
cate that the approximation of unstable subspace via singular value decomposition 
is preferable to that of Schur decomposition, the theory at present holds only for 
Schur decomposition. Future work should concentrate on the mathematical analysis 
needed to formulate similar ideas for singular value decomposition; a clearer math- 
ematical understanding should enable the development of an efficient strategy for 
systematic detection of periodic orbits in high- dimensional systems. 

We construct STs by approximating the local stretching rates of the system at 
the initial seed x. As we evolve the associated flow however, we can wander far away 
from the initial seed, in that case the ST matrix is no longer valid due to its local 
nature. One might try to improve the applicability of the method by constructing 
transformations adaptively, which amounts to continually updating our approxima- 
tion of the unstable subspace. However, since we approximate the unstable subspace 
at X by evaluating the singular value decomposition at its preimage, i.e. for 
some time t, we must be able to integrate backwards in time. Unfortunately, for 
dissipative systems such as the KSE, evolving the solution for negative times is not 
possible. However, this is not the case for reversible systems, for example Hamilto- 
nian systems, in that case we can integrate either forward or backwards in time in 
a straightforward manner. 

Another possible avenue of exploration would be the application of ST approach 
to Hamiltonian systems, which are most relevant for the study of quantum dynam- 
ics of classically chaotic systems and celestial mechanics. This should allow for 
further improvements - particularly theoretically - to be achieved through the con- 
sideration of the symmetry properties of the Hamiltonian systems. Indeed, since 
the choice of the transformations is directly related to the eigenvalues of the orbit 
monodromy matrix, one of the obvious properties that can be exploited in this case 
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is the time-reversible nature of the Hamiltonian systems, which manifests itself in 
the symmetry of the eigenvalues of the corresponding monodromy matrix. Other 
symmetry properties related to the symplectic structure could also be explored. 

An important consideration when applying our method to Hamiltonian problems 
is that for a given energy the search for UPOs take place on the corresponding 
energy surface. However, in the determination of UPOs we solve an associated 
flow which, in general, will cause us to move away from the energy surface. One 
possible approach to avoid this, would be to project out any part of the vector 
field transverse to the surface and to solve the resulting system of ODEs, thus 
remaining on the energy surface. Alternatively, we could add an extra equation, 
which would act so as to force the associated flow onto the energy surface in an 
asymptotic manner. This could work, for example, in a similar fashion to penalty 
merit functions in nonlinear programming problems, where one typically has to 
make a tradeoff between minimising the function of interest and remaining on the 
manifold of feasible points. 



Appendix A 



Frequently used notation 



[p\ Period of a discrete system 6 

[/] Discrete map 7 

[g] F{x) -X 7 

[n] System dimension 7 

\U] A discrete dynamical system 7 

[S] Associated flow 7 

[C] Stabilising transformation 7 

[Csd] Stabilising transformations proposed by Schmelcher and Diakonos 7 

[G] The Jacobian of ^'(x) 8 

n X n identity matrix 24 

[DfP] The jacobian of 24 

[Dg] Alternative notation for the Jacobian oi g{x) 25 
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[T] Matrix transpose 27 

[0*] The flow map of a differential equation 30 

[T] Period of a continuous system 30 

[J] The Jacobian for the flow J = D(f)^{x) 30 



Appendix B 



Abbreviations 



[UPO] Unstable periodic orbit 1 

[SD] Schmelcher-Diakonos 6 

[ODE] Ordinary differential equation 8 

[PDE] Partial differential equation 9 

[ST] Stabilising transformation 15 

[PSS] Poincare surface of section 15 

[BW] Biham-Wenzel 20 

[CB] Characteristic bisection 23 

[NR] Newton-Raphson 24 

[GN] Gauss-Newton 34 

[LM] Levenberg-Marquardt 35 

[CHM] Coupled Henon map 55 
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[SchD] Schur decomposition 62 

[SVD] Singular value decomposition 63 

[KSE] Kuramoto-Sivashinsky equation 68 

[FFT] Fast Fourier transform 69 



Appendix C 



Derivation of the kicked double 

rotor map 



The double rotor map is a four-dimensional discrete system, physically it describes 
the dynamics of a double rotor under the influence of a periodic kick; see Figure 
C.l. Note that our discussion closely follows that given in [32], we start by giving 
a description of the mechanical device that is the double rotor, before deriving the 
equations of motion and the corresponding map. 

The double rotor is made up of two thin, massless rods, connected as in figure 
C.l. The first rod, of length /i, pivots about the fixed point Pi, whilst the second 
rod, of length I2, pivots about the moving point P2- The position at time t of the 
two rods is given by the angular variables 0i{t) and 62{t) respectively. A mass mi is 
attached to the first rod at P2, and masses 7712/2 are attached at either end of the 
second rod. Friction at Pi is assumed to slow the first rod at a rate proportional to 
61, and friction at P2 slows the second rod at a rate proportional to 62 — 9i. The end 
of the second rod, marked K, receives a periodic kick at times t = T, 2T, 3T . . . , the 
force of which is constant. 

We may construct the equations of motion for the system by noting that the 
kick represents a potential energy of the system which is given by 



here f(t) consists of periodic delta function kicks, i.e. f(t) = J2k /o^(^ ~ kT) . The 
kinetic energy of the system is easily seen to be 



V{t) = {h cosOi + k cos 62) f{t) 



(C.l) 



m = -{hel + i2el) 



(C.2) 
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Figure C.l: The double rotor under the influence of a periodic kick. 

where Ii = (mi + 7712)/^, I2 = ^2/2 and 6i = d9i/dt for i = 1,2. Using the Lagrangian 
formulation, L = K — V, with 

l(dL\_dL^_dF 

dt \dqj % de, ' ^ ■ ' 

where Rayleighs dissipation function is given by F = ^uiIiOi + |z/2-^2(^2 — ^1)^, 
yields the following system of ODEs 

Ii^ = Ii/(t)sin^i-z/ili^i + z/2l2(^2-^i), (C.4) 

l2-^ = l2/(t)sin^2-^^2l2(^2-^^l), (C.5) 

at 

here ui, 1/2 are the coefficients of friction. 

Since the kicks are treated as being instantaneous, i.e. f(t) = 0, except for times, 
t, which are integer multiples of the period, the equations of motion (C.4) - (C.5) 
reduce to the following linear system for t G {{n — 1)T, nT), n G Z+, 



d^i/dt \ _ f -(z/i + z/2) z/2 \ / ^1 
d^2/dt I \ z/2 -1^2 



(C.6) 



Given the initial conditions (6*1(0), ^2(0))""", Eq. (C.6) admits the following solution 

«2(i) / \ «2(0) ' 
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where 



UuUi2e + M2lM22e ''^^ 



UiiUue 



— Sit 



-S2t 



-Sit 



1X22^ 



-S2t 



(C. 



Si, S2 are eigenvalues of the matrix in Eq. (C.6) with corresponding orthonor- 
mal eigenvectors ui = (^11,^12)"'", U2 = (^21, M22)"'"- For the initial condition 
(^i(O), 6'2(0)), we may integrate Eq. (C.7) to determine the position of the rods 




M{t) 




(C.9) 



here M{t) = J^L{^)d^. Note that Eqs. (C.8) - (C.9) completely describe the dy- 
namics of the double rotor for t G {{n — 1)T, nT), n G Z"*". 

At t = T the angular velocity of each rod changes instantaneously, thus the 
angular velocity 9i of each rod is discontinuous. Denote the limits from the left and 
right as 6i{T~), 6i{T~^) respectively, then the discontinuity is given by 



9,{T+)-9,{T- 



^sine,(T), 



1,2. 



(C.IO) 



The position however, will vary continuously for all time, i.e. 6i{T^) = 6i{T~), 
i = 1,2. In this way the solution of (C.4) - (C.5) is a composition of the solution of 
the linear system (C.6) and the periodic kicks at t = T, 2T, .... Thus to understand 
the dynamics of the double rotor it suffices to study the following map obtained 
from Eqs. (C.7) - (C.IO), 



^(n+1) 



9, 



5(n+l) 



M{t) 
L{t) 



( \ 

V ^2 / 



+ 



in) 



\ 



hfo 



(n+1) 



^ sin e\ 
^sin^5"+') 



(C.ll) 



which give the position and the angular velocity of the rods immediately after each 
kick. 



Appendix D 



A MODICUM OF LINEAR ALGEBRA 



Below we define some of tlie tlieory from linear algebra used throughout the current 
piece of work. Proofs are omitted for brevity but can be found in the books on matrix 
analysis by Horn and Johnson from which the results have been taken [45, 46]. In 
what follows we denote by M„ the set of n x n, possibly complex valued matrices. 
Further we denote by A* the conjugate transpose of A, this of course equals the 
usual transpose for real matrices. 

We begin with two theorems concerning the matrix decompositions which have 
been used so frequently in our discussions 

Theorem D.l. If G E M„, then it may be written in the form 

G = QB (D.l) 

where B is positive semi-definite and Q is unitary. The matrix B is always uniquely 
defined as B = [GG*)^^"^ ; if G is nonsingular, then Q is uniquely determined as 
Q = B^^G. If G is real, then B and Q may be taken to be real. 

The factorisation D.l is known as the polar decomposition; an important point to 
note here is the uniqueness of the two factors in the case that G is nonsingular. 
A perhaps even more useful decomposition is that of singular value decomposition 
(SVD). Its applications are far reaching and include for example, the study of linear 
inverse problems, uses in signal processing and many areas of statistics. 

Theorem D.2. If G E Mp^n has rank k, then it may be written in the form 



G = VSW* 



(D.2) 
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where V and W are unitary. The matrix S = [sij] G Mp „ has Sij = for all i ^ j , 
and Sii > S22 > ■ ■ • > > Sk+i^k+i = ■ ■ ■ = Sgg = 0, where q = min(p, n). The 
numbers {sa} = {si} are the nonnegative square roots of the eigenvalues ofGG*, and 
hence are uniquely determined. The columns of V are eigenvectors of GG* and the 
columns ofW are eigenvectors of G*G. Ifp < n and if GG* has distinct eigenvalues, 
then V is determined up to a right diagonal factor D = diag(e*^i, . . . , e*^") with all 
9i e M; that is, if G = ViSW^ = ViSW^, then V2 = ViD. Ifp < n, then W is never 
uniquely determined; if p = n = k and V is given, then W is uniquely determined. 
If n <p, the uniqueness ofV and W is determined by considering G* . if G is real, 
then V , S and W may all be taken to be real. 

The two decompositions admit the following relation 

G = QB = VSW* (D.3) 
= VW*WSW* (D.4) 

thus we have that Q = VW* and B = WSW*. In constructing the stabilising 
transformations of Chapter 3, we use the SVD of the matrix G, rather than its 
polar decomposition, in our numerical calculations. 

Below we will give the details of the theorems of Lyapunov and Sylvester which 
where used in the proof of Corollary 3.1. However, before doing so we find it instruc- 
tive to present a couple of useful definitions of what we believe to be nonstandard 
linear algebra. The first notion is that of the inertia of a general matrix in M„ 

Definition D.l. If A e M„, define: 

i+{A) = the number of eigenvalues of A, counting multiplicities, with positive real 
part; 

i-{A) = the number of eigenvalues of A, counting multiplicities, with negative real 
part; and 

io{A) = the number of eigenvalues of A, counting multiplicities, with zero real part. 
Then, i+{A) + i-{A) + io{A) = n and the row vector 

i{A)^[i^{A),i4A),i){A)] 



is called the inertia of A. 
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This leads us nicely onto our second definition, which simply translates our idea of 
stability in terms of a dynamical system, into the language of linear algebra 

Definition D.2. A matrix A e Mn is said to be positive stable if i{A) = [n, 0,0], 
that is, i+{G) = n. 

The matrices which we study in this work are typically derived from some dynamical 
system of interest. In this case, stability is defined in terms of matrices who possess 
only eigenvalues with negative real parts, however, it is convention in linear algebra 
to discuss positive matrices, and we shall stick to convention. It should be clear 
that to determine that a matrix A is stable, it suffices to show that the matrix —A 
is positive stable. Our last definition of this section concerns the idea of congruence 
of two matrices 

Definition D.3. Let A,B E Mn he given. If there exists a nonsingular matrix S 
such that 

A = SBS* (D.5) 
then A is said to be *congruent ("star- congruent") to B. 

The defining property for the equivalence class of *congruent matrices introduced 
above is that their inertia remains constant. That this is so is given by the following 
theorem which is known as Sylvester's law of inertia 

Theorem D.3. Let A,Be Mn be Hermitian matrices. There is a nonsingular 
matrix S G M„ such that A = SBS* if and only if A and B have the same inertia, 
that is, the same number of positive, negative, and zero eigenvalues. 

We finish the section with the following result due to Lyapunov 

Theorem D.4. Let A g M„ be given. Then A is positive stable if and only if there 
exists a positive definite G G Mn such that 

OA + A*G = H (D.6) 

is positive definite. Furthermore, suppose there are Hermitian matrices G,H E Mn 
that satisfy (D.6), and suppose H is positive definite; then A is positive stable if and 
only if G is positive definite. 

The theorem gives a nice relation between the class of positive definite matrices and 
the class of positive stable matrices. 



Appendix E 



Exponential time differencing 



Stiff systems of ODEs arise naturally wlien solving PDEs by spectral methods, and 
their numerical solutions require special treatment if accurate solutions are to be 
obtained efficiently. In this appendix we describe one such class of solvers, the 
Exponential time differencing (ETD) schemes, they where first used in the field 
of computational electrodynamics, where the problem of computing electric and 
magnetic fields in a box typically result in a stiff system of ODEs; see [12] and 
references therein for further details. In this discussion we concentrate on the Runge- 
Kutta version of these schemes, in particular, we look at a modification of the ETD 
fourth-order Runge-Kutta method presented by Kassam and Trefethen [51]. 
Let us represent our PDE in the following form 



here C and A/" are linear and nonlinear operators respectively. Applying a spatial 
discretisation yields the following system of ODEs 



To derive the ETD methods, we begin by multiplying Eq. (E.2) by the integrating 
factor e""""*, and integrating over one time step from t = t„ to t = tn+i = tn + h to 
obtain 



^0 

this formula is exact, and the different ETD methods result from the particular 
choice of approximation for the integral in Eq. (E.3). For example, if we assume 



Ut = Cu + A/'(m, t). 



(E.l) 



Ut = Lu + N(m, t). 



(E.2) 




(E.3) 
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that N is constant, i.e. N = N(u„) + 0{h), over a single time step then we obtain 
the ETDl method 

Un+i = u^e^^ + L-i(e^^ - l)N(n„), (E.4) 

which has a local truncation error /i^N/2. In [12] Cox and Matthews present a host 
of recurrence formulae that provide higher-order methods, as well as introducing 
the set of methods based on Runge-Kutta time-stepping which they name ETDRK 
schemes. 

In our work we have used the fourth order scheme, known as ETDRK4, the 
derivation of which is nonstandard according to Cox and Matthews. For this reason 
we simply present the formulae and refer the interested reader to [12] 

an = e^^/2M„ + L-i(e^'^/2-I)NK,t„), (E.5) 
hn = e^'''^Un + 'L~\e^'''^ - m{an.tn + h/2), (E.6) 
c„ = e^^/V + L-^(e^'^/2-I)(2N(6„,t„ + V2)-NK,t„)), (E.7) 
M„ = e^^/2 + /i-2L--^{[-4-L/i + e^'^(4-3L/i+(L/i)2)]N(M„,tn) (E.8) 
+2[2 + L/i + e^^(-2 + L/i)] (N(a„, t„ + h/2) + N(6„, t„ + h/2)) 
+ [-4 - 3L/i - (L/i)2 + e^^(4 - L/i)]N(c„, tn + h)]. 

Unfortunately, the stated method is prone to numerical instability, the source of 
which can be understood by examining the following expression 

f{z) = (E.9) 

The accurate computation of this function is a well known problem in numerical 
analysis, and is further discussed and referenced in the paper by Kassam and Tre- 
fethen; the difficulties are born from cancelation errors which arise for small z. 

To understand why Eq. (E.9) relates to the ETDRK4 method it is useful to 
examine the coefficients in the square brackets of the update formula 

a = h-^L-^[-4-Lh + e^''{4-3Lh + (Lh)% 

(3 = h-'^L-'^[2 + Lh + e^''{-2 + Lh)], (E.IO) 

7 = h-^L-^[-4-3Lh- {Lhf + e^''{4-Lh)]. 

These coefficients are actually higher-order variants of Eq. (E.9), and thus, are sus- 
ceptible to cancelation errors. Indeed, all three terms suffer disastrous cancelation 
errors when the matrix L has eigenvalues close to zero. Cox and Matthews knew 
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of this problem in their work, and proposed a cutoff point for small eigenvalues, 
whereafter the coefficients would be computed by Taylor series. However, two prob- 
lems arise: firstly, there may exist a region where neither the formula nor the Taylor 
series are accurate, and secondly, it is not obvious how to extend these ideas in the 
case that L is not diagonal. 

In order to obtain numerical stability the accurate computation of the coefficients 
(E.ll) is of paramount importance. With this in mind Kassam and Trefethen suggest 
evaluating the coefficients via an integral over a contour in the complex plane. Recall 
from Cauchy's integral formula that a function f{z) may be computed as follows 



where F is a contour that encloses z and is separated from 0. The application of 
Cauchy's formula in the case that z is a matrix is straightforward 



notice that the term (t — z)~^ has been replaced by the resolvent matrix, and that 
we now choose F so that it encloses all the eigenvalues of L. 

When the above stabilisation procedure is applied to the ETDRK4 method the 
result is a stiff solver which works equally well, whether the linear part of Eq. (E.2) is 
diagonal or not, it is extremely fast and accurate, and allows one to take large time- 
steps. For example, in the computation of the finite-time Lyapunov exponents of the 
KSE in Chapter 4, we where able to use step-size h = 0.25, the same computation 
using the eighth-order RK method due to Dormand and Prince [37] needed to use 
h = 2x 10"^ 




(E.ll) 




(E.12) 



Appendix F 



Numerical calculation of 
Lyapunov spectra 



In this section we briefly remind ourselves of the idea of Lyapunov exponents, before 
discussing an algorithm first introduced by Benettin et al [2] in order to efficiently 
determine them; in Chapter 4 we use this algorithm to determine the largest Lya- 
punov exponents for the KSE. 

The concept of Lyapunov exponents is an important notion for dynamical sys- 
tems, particularly in applications, where it is often used as a criteria for determining 
the existence of chaos. In the following we give a description in the case of a dis- 
crete dynamical system, x^+i = f{xk)- In that case the Lyapunov exponents give a 
description of the average behaviour of the derivative, Df, along the orbit of some 
initial point xq, i.e. {xq, /(xq), P{x), . . .}. 

Write Eq for the unit ball in ra-dimensional phase space centred at xq, and define 
successive iterates by 

Ek+i = Df{xk)Ek. (F.l) 

Note that each E^ is an ellipsoid. Now, let aj^k denote the length of the jth largest 
axis of Eki then we define the jth Lyapunov number of / at xq to be 

Lj= lim(a,•fc)l/^ (F.2) 

when the limit exists. The Lyapunov exponents are then given by the natural loga- 
rithms 

Xj = \nLj. (F.3) 
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The trajectory of the point Xq is called a chaotic trajectory if (1) the trajectory is 
bounded and is not asymptotic to an equilibrium position of /, and (2) / has at 
least one positive Lyapunov exponent. 

Computation of Lyapunov exponents 

As noted an efficient algorithm for determining the Lyapunov exponents of chaotic 
orbits has been put forward by Benettin et al and is as follows: starting with an 
orthogonal set of unit vectors . . . , Define 

Vi = Df{x)ui 

where i — 1, . . . ,n. Applying the Gram-Schmidt algorithm we compute a set of 
orthogonal vectors, Wi {i — 1, . . . ,n), such that wi — vi, and for i > 1, the Wi are 
defined inductively as follows 

i-l 

Wi = Vi-^< Vi, Uj > Uj, (F.4) 

here Uj = Wj/\\wj\\, where || ■ || denotes the Euclidean norm, and < •, • > denotes 
the corresponding inner product. 

The new set of vectors, uj, are approximations of the directions of the ith axis 
of Ek+ii further, the approximate ratio of the length of the ith axis of Ek+i to that 
of Ek is given by = \\wi\\. To determine our approximation this operation is 
repeated on each iterate of the map /. Writing rj j for the ratio at the jth iterate 
of / we obtain the following approximation to the ith Lyapunov exponent 

k 

The approximation improves with increasing k. 
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